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ABSTRACT 

In  the  present  lecture  some  basic  problems  of  state  and  solving  of  radiation  heat  transfer  equation  as 
applied  to  radiation  modelling  in  shock-tubes  and  entry  flow  are  discussed.  The  lecture  contains  five 
parts. 

In  the  first  part  the  radiation  heat  transfer  equation  and  the  general  definitions  of  the  radiation  heat 
transfer  theory  are  presented.  The  definitions  introduced  in  the  first  part  are  widely  used  in  other  parts. 

The  second  part  of  the  lecture  ( Chapter  3)  presents  four  methods  which  are  used  for  solution  of  radiation 
heat  transfer  problems  in  different  aerospace  applications. 

Well  known  problem  of  radiation  heat  transfer  in  heated  gases  and  low -temperature  plasma  in  view  of 
atomic  and  molecular  rotational  lines,  as  well  as  of  the  vibrational  molecular  bands,  are  discussed  in  the 
third  part.  Some  untraditional  models  of  radiation  heat  transfer  are  presented  in  this  part. 

The  fourth  part  of  the  lecture  is  dedicated  to  the  Monte-Carlo  algorithms  which  are  also  in  common  use  in 
aerospace  problems,  especially  at  prediction  of  emissivities  of  heated  and  radiating  volumes  of  light¬ 
scattering  gases.  Several  traditional  and  novel  algorithms  are  presented  in  this  part. 

The  final  part  ( Chapter  6)  of  the  lecture  presents  examples  of  using  of  numerical  simulation  methods 
described  in  previous  parts  at  solution  of  various  radiation  heat  transfer  problems  in  aerospace 
applications.  It  should  be  stressed  that  the  methods  of  the  radiation  heat  transfer  integration  and  various 
numerical  simulation  results  presented  in  the  lecture  were  used  and  obtained  by  the  author  at  solution  of 
concrete  problems  of  aerospace  physical  gas  dynamics. 
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1.0  INTRODUCTION 

The  term  “thermal  radiation”  (“heat  radiation”)  means  electromagnetic  radiation  of  atomic  and  molecular, 
as  opposed  to  nuclear,  origin.  Such  radiation  is  emitted  by  matter  in  a  state  of  thermal  excitation,  thus 
accounting  for  the  designation  of  the  radiation  as  thermal.  The  energy  density  of  this  type  of  radiation  is 
given  by  the  Planck  formula  for  black  body  radiation.  More  generally  the  radiation  energy  distribution  is 
described  by  a  kinetic  or  transport  equation,  referred  to  historically  as  the  equation  of  radiative  transfer. 

The  importance  of  thermal  radiation  in  physical  problems,  and  particularly  in  problems  of  physical  gas 
dynamics  in  different  aerospace  applications,  increases  as  the  temperatures  is  raised.  At  moderate 
temperatures  (104  <  T  <105K),  the  role  of  radiation  is  primary  one  of  transporting  energy  in  gases  and 

~  ~  c  /r 

plasmas  by  radiative  processes.  At  higher  temperatures  (  T  >10  -MO  K),  the  energy  and  momentum 
densities  of  radiation  field  may  become  comparable  or  dominate  the  corresponding  fluid  quantities.  As  a 
rule,  hydrodynamics  with  explicit  account  of  the  radiation  energy  and  momentum  contributions  constitutes 
subject  of  investigation  of  “high-temperature  radiation  hydrodynamics”.  More  general  definition  of  radiation 
gasdynamic  implies  consideration  of  coupled  radiative  and  gasdynamic  processes.  In  the  partial  case  of  the 
weak  radiation-gasdynamic  interaction  there  is  possibility  to  study  radiative  and  gasdynamic  phenomena 
separately.  For  example,  such  a  case  is  realized  for  spacecraft  at  entry  velocities  up  to  6-^8  km/s. 

Radiative  gas  dynamics  (RadGD)  is  the  directions  of  physical  gas  dynamics  which  is  connected  to  such 
challenging  sciences  as:  astrophysics,  physics  of  stars  and  Sun,  research  of  a  structure  of  substance 
(atomic  and  molecular  spectroscopy),  interaction  of  laser  radiation  and  high-energy  beams  with  materials, 
plasma  generators,  rocket  engines  (of  chemical,  plasma,  electric,  nuclear  or  laser  types),  spacecraft's 
thermal  protection,  heat  exchange  in  steam  boiler,  in  aircraft  and  rocket  engines,  in  working  volumes  of 
various  power  installations  (including  nuclear). 

Figure  1.1  shows  hierarchical  division  of  the  Radiative  gas  dynamics.  In  order  to  solve  any  RadGD 
problem  there  is  a  necessity  to  create  a  radiative  model  of  a  gas  or  plasma.  The  radiative  model  is  defined 
as  the  set  of  optical  model  and  radiation  transfer  model. 

The  optical  model  includes  spectral,  group  and  integral  absorption,  emission  and  scattering  coefficients, 
which  are  in  turn  based  on  cross-sections  (or  probabilities)  of  elementary  radiative  processes  predicted  by 
quantum  mechanics  and  quantum  chemistry.  The  absorption,  emission  and  scattering  coefficients  can  be 
determined  only  with  use  data  on  distribution  functions  for  atomic  and  molecular  particles,  and  also  on  their 
energy  states.  Thermodynamics  and  statistical  physics  provide  all  necessary  information  for  these  purposes. 

The  radiation  transfer  model,  composing  the  second  part  of  the  radiative  model,  is  based  on  the 
thermodynamics  and  statistical  physics,  and  is  designed  for  prediction  such  characteristics  of  a  radiation 
field  as  the  radiation  energy  density  U,  the  radiation  flux  W  ,  and  the  divergence  of  the  radiation  flux 
Qrad  =  di vW  .  A  spectral  region  of  the  electromagnetic  radiation,  which  is  of  practical  interest  for  various 
aerospace  applications  ranges  from  □  0.05  pm  until  □  20  pm.  This  spectral  region  is  divided  on  the 
following  sub  regions: 

•  AX  «  0.05  -r-  0.4  pm  is  the  ultraviolet  region; 

•  AX  «  0.4  -r-  0.7  pm  is  the  visible  region; 

•  AX&  0.7  -r-  20.0  pm  is  the  infrared  region, 

where  AX  is  the  wavelength  interval.  A  radiation  emitted  by  matter  in  this  spectral  region  is  called  as  heat 
radiation. 
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Figure  1.1 :  Basic  scheme  of  Radiative  Gas  Dynamics 


2.0  RADIATION  HEAT  TRANSFER  EQUATION  AND  GENERAL 
CHARACTERISTICS  OF  HEAT  RADIATION  TRANSFER 

2.1  Classification  of  spectral  optical  models 

Classification  of  the  spectral  optical  models  was  introduced  in  [1]: 

1)  Optical  model  /^(T^cp)  is  spectral,  quasi-spectral,  multigroup  and  combined  models  of  absorption 
coefficients  representing  in  the  table,  analytical,  graphical  or  in  a  computer  module  form.  Here  co  is 
the  Wavenumber;  T  is  the  temperature;  cp  is  set  of  the  parameters  (pressure  or  density,  concentration 
of  the  chemical  components).  In  the  case  of  a  nonequilibrium  medium  instead  of  T  we  should  use  set  of 
effective  temperatures  (electronic  temperature,  vibrational  temperature,  rotational  temperature). 

2)  Spectral  ( line-by-line )  model  (T,cp)  is  the  function  with  continuous  and  lines  structure  without 
any  smoothing. 

3)  Multigroup  model  is  the  smoothing  spectral  model  in  a  set  of  spectral  ranges  Acoz .  We  suppose  in 
the  each  spectral  range  Aco,  the  absorption  coefficient  is  independent  of  wavenumber. 

4)  Quasi-spectral  ( quasi-continuum )  model  is  the  smeared  rotational  line  model  for  a  set  of  spectral 
ranges.  The  value  Aco,  in  this  model  must  be  larger  than  maximal  values  of  the  rotational  lines  widths 
(that  is  not  less  than  25  -f  50  cm-1).  As  a  rule,  quasi-spectral  model  is  included  to  multigroup  model. 
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5)  Total  absorption  model  is  the  Planck  mean  absorption  coefficient,  or  (and)  the  Rosseland  mean 
absorption  coefficient,  or  (and)  the  Chandrasekhar  mean  absorption  coefficient. 

6)  Combined  model  is  the  sum  of  the  spectral  line  absorption  and  the  multigroup  absorption  model  (as 
a  rule  the  selective  atomic  lines  absorption  on  the  continuum  or  quasi-continuum  background). 

7)  Radiative  heat  transfer  model  is  the  set  of: 

•  conditions  (the  thermodynamic  conditions  in  a  medium,  radiation  heat  transfer  boundary 
conditions,  the  spectral  resolution  of  numerical  calculations); 

•  equation  of  radiation  heat  transfer; 

•  mathematical  method  using  for  solving  this  equation. 

8)  Radiative  model  is  the  set  of  the  optical  model  and  the  radiation  heat  transfer  model. 

9)  Optimum  radiative  model  is  the  radiative  model  containing  minimum  number  of  the  spectral  groups 
and  in  the  same  time,  permitting  to  obtain  physical  adequately  results. 

Figure  2.1  shows  example  of  air  low-temperature  plasma. 


Wavenumber,  1/cm 

Figure  2.1 :  Line-by-line  and  group  spectral  model  of  absorption  coefficient  of  air  at  7=  10000  K  and  p  =  1  atm 

2.2  General  notations  of  the  radiation  heat  transfer  theory 

The  steady-state  radiation  heat  transfer  equation,  which  is  taking  into  consideration  scattering  processes 
has  the  following  form  [2,  3]: 

—  ^ +  KV7V  (  j,Q)  +  av/v  (s,Q)  =  jv  (s)  +  ^  J  pv(s,Cl'  ^Cl)jv(Cl')dCl'  ,  (1) 

1  4  71 

where  is  the  spectral  radiation  intensity;  kv(s)  is  the  spectral  absorption  coefficient;  ov(s)  is 

the  spectral  scattering  coefficient;  jv(s)  is  the  spectral  emission  coefficient;  is  the  spectral 
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phase  function  for  scattering  (the  scattering  indicatrix),  s  is  the  physical  co-ordinate  along  spatial 
direction  Q  .  Angular  variables  for  unit  vector  Q  are  presented  in  Figure  2.2. 


Figure  2.2:  Angular  co-ordinates  for  integration  of  a  spectral  radiation  intensity  on  solid  angle 

The  following  type  of  the  equation  is  especially  often  uses: 


(2) 


where 


is  the  Planck  intensity  (the  spectral  intensity  of  the  black  body  emissivity);  h,k  are  the  Planck  and 
Boltzman  constants;  c  is  the  speed  of  light. 

Equation  (2)  is  correct  for  the  Local  Thermodynamic  Equilibrium  (LTE),  and  does  not  take  into  account 
scattering  processes. 

The  following  units  for  the  electromagnetic  spectra  measurement  are  in  common  use:  X  is  the  wavelength 
(in  microns,  Angstroms,  nanometers;  1  pm  =  1CT6  cm,  1A  =  1CT8  cm,  1  nm=  10-9m);  v  is  the  frequency  (in 
s-1);  co  is  the  wavenumber  (in  cm-1).  These  units  are  connected  with  the  following  relations  (for  radiative 
processes  in  vacuum):  v  =  cco,  c  =  vX.  This  allows  to  write  the  following  correlations  for  the  Planck 
function:  Zvdv  =  /^dco,  /Vdv  =  ~J\dk .  To  recalculate  the  wavelength  to  the  wavenumber  one  can  use  the 
following  relation:  =104/cocm-i  . 

The  wavelength  and  wavenumber  dependencies  of  the  Planck  function  at  different  temperatures  are  shown 
in  Figures  2.3.  These  data  together  with  the  Wien  displacement  law  XmaxT  =  2897.8  pm-K  are  very 
suitable  for  express  analysis  of  the  radiation  heat  transfer  problems.  Here  Xmax  is  the  wavelength,  in  pm, 
at  which  the  emissive  power  /  is  a  maximum  for  a  given  temperature  T  (in  K). 

The  spectral  radiation  intensity  allows  define  the  following  general  spectral  and  integral  characteristics  of 
the  radiation  field: 

C/qj  =  —  J  J&  ^,Q^dQ  is  the  spectral  energy  density; 


4tt 
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4r,„,co=(^V”)  =  f  4,(s,^)p-«)dn  is  the  spectral  hemispherical  flux; 

471 

00  00  00 

u  =  \  Uda  ,W  =  J  Wd© ,  qRn  =  J  qR  n  (i) dco , 

0  0  0 

are  the  integral  energy  density,  integral  radiation  flux,  and  integral  hemispherical  flux. 


The  Planck  function,  W/(ster*cm**2*mcm) 


Figure  2.3:  The  wavelength  (left)  and  wavenumber  (right)  dependence  of  the  Planck  function  at  different 

temperatures 


The  radiative  energy  balance  in  a  gas  and  plasma  can  be  calculated  as  follows: 


J  QV/dQ  =  J  K/M[r(s)]dn-  J  K/m(5,Q)dQ,  (3) 

471  471  471 

or  V-W(0=47iK(B7/7(0-Kj7m(5,Q)dQ, 

471 

where  the  first  summand  in  the  right  hand  side  corresponds  to  the  spectral  emission,  and  the  second  one 
corresponds  to  the  spectral  absorption. 

If  we  integrate  Equation  (3)  by  spectrum,  then 


( V  •  W^j  =  47idco  -  c  J  Kf/dco ,  (4) 

o 

00 

or  (V  •  W)  =  4k pdT4  - c J  kw  U0) dco  ,  W/cm3, 

0 


(*  CJ  —12  W 

where  /dco  =  —T,  a  =  5.67  x  10  ,  — - — -  is  the  Stefan-Boltzmann  constant; 

J  n  Cm2K4 
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kp  =  |  ic/dcoy/J  /doo  is  the  Planck  mean  absorption  coefficient. 
Other  significant  characteristics  of  radiation  heat  transfer  are: 


1)  The  group  absorption  coefficients  are  introduces  as  following: 

1  j  ^Acog 

k  =  kAco  = -  [  Kdco  = -  V  kAco^  is  the  group  absorption  coefficient; 

*  Aco  Acg  ri 

Am  I—1 


(5) 


Wa, 


k P  g  =— —  ^  k/Aco(  is  the  group  Planck  mean  absorption  coefficients; 

i=l 


(6) 


Ns 

kp  =  > 

g= i 

where  NAm  is  the  number  of  the  spectral  points  inside  g-th  spectral  group,  Ng  is  the  number  of  the  group, 
Am,  is  the  elementary  spectral  region.  Actually,  Equations  (5)  and  (6)  introduce  the  line-by-line  model. 

2)  A  photon  mean  free  path 


Lr  = 


ou 

=  -5-f 

A-rp?)  J 


1  8J, 


b,  co 


4dTiJnKm  dT 


0  © 
N Am  „ 


4  dT 


n  y?  1 


dco,cm, 


( A 

OJb,  co 


(8) 


i=l  ^Aco, 


sr 


Af 


Aco(- ,  LR-^LRg. 


A 


?=l 


3)  The  Rosseland  mean  absorption  coefficients  are  calculated  as  following: 


4) 


1 


1 


Kd  —  - 


-R,g 

The  radiative  thermal  conductivity  is  defined  as 

Qr  =  -kfiVT,  TR  =  —  dT . 

3  Kp 


(9) 


L 

This  approximation  is  correct  only  for  x^  □  1 ,  where  =  J  K^ds  is  the  optical  depth;  L  is  the  length  of 

o 


the  optical  path. 


5) 


The  spectral,  group,  and  integral  emissivities  can  be  defined  in  the  following  form  only  for  the  LTE 
approximation: 


Qemy  4t11C/ by 


w 


cm3  •  jlx 


2Am  =4KP,gar>  d 


cem  =4k6T\^j 

cm 


(10) 


Because  the  Planck  function  is  isotropic  one,  the  integrated  flux  of  the  absolutely  black  body  can  be 
calculated  as  follows: 
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%,n  (T)  =  \'{)  do)j,^ d(pj,f2 4„SinOCosOd0  =  dT4 . 


(11) 


For  estimation  and  comparison  of  emissivity  of  various  mixes  of  high-temperature  gases  a  spherical 
emissivity  factor  is  used 


(12) 


In  particular,  the  emissivity  factor  of  a  half-spherical  homogeneous  volume  of  a  radius  R  can  be  calculated 


as 


R  2n  n/2  , 

ev,sph  =(^»r1 2fdrj  d(P  J  Ka)4a)exP(-Ka)r)Cos0Sin0d0=  ^T^  =  1-exP(-Kco/?)'  (l3) 

0  0  0 

The  emissivity  factor  of  a  flat  homogeneous  layer  is  determined  by  the  following  expression 


L  2k  jc/2 

W  =  (^/xo )  ' jcu  j  d(P  {  K aJba exp 
0  0  0 


K...X  A 


CosB 


Sin0d0  = 


TlJ, 


■  1-2E3(k,L) 


(14) 


bco 


1 

where  En  (x)  =  J exp(-x/p)pw_2dp  is  the  integro-exponential  function  of  the  n- th  order,  L  is  the  depth  of 
o 

the  flat  layer. 


An  integral  emissivity  factors  are  calculated  as  following: 


eAco  =  ~A  J  £»dt0’  £  =  if £»dc0  • 

Aco  0 


(15) 


There  is  relation  between  spectral  absorption  coefficient  and  spectral  emission  coefficient  in  the 
LTE  conditions,  which  is  based  on  the  Kirchhoff  law:  /w  =  k . 


3.0  METHODS  OF  INTEGRATION  OF  RADIATION  HEAT  TRANSFER 
EQUATIONS 

3.1  The  Pn  (Spherical  Harmonic)  approximation  [2,  4] 

A  set  of  approximations  to  the  equation  of  radiation  transfer,  which  is  capable  of  estimating  the  solution  to 
the  equation  of  transfer  to  within  any  desired  accuracy  criteria,  is  considered  in  the  chapter.  This  is  the  so- 
called  PN-approximation  of  the  spherical  harmonic  method,  with  N  denoting  the  order  of  the 
approximation.  For  N  infinite,  the  PN-approximation  is  the  exact  solution  to  the  radiation  transfer  equation. 


The  PN-approximation  will  be  considered  for  two  calculation  cases,  which  are  of  practical  interest  for 
aerospace  applications. 

1)  Radiation  heat  transfer  in  infinite  plane  inhomogeneous  layer  with  light-scattering  media  (ID 
radiation  heat  transfer  problem). 

2)  Radiation  heat  transfer  in  axially  symmetric  two-dimensional  volume  (2D  radiation  heat  transfer 
problem). 
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3.1.1  Plane-parallel  calculation  case 

Let  us  consider  radiation  heat  transfer  equation  for  the  ID  calculation  case  shown  in  Figure  3.1,  which  is 
described  by  the  following  equation: 

~  +  /(t^)  =  (1_5>)4(t)+^  J  I-041'  (i) 

with  the  boundary  conditions 

x  =  0,  n>0:  /(0,n)  =  /0+(n),  (2) 

t  =  t#,  n<0:  y(x//,n)  =  y^(n).  (3) 

where  /(x,|u)  is  the  spectral  intensity  of  heat  radiation  (index  of  spectral  dependence  is  omitted);  x  is  the 
spectral  optical  depth  dx  =  Kcbc;  k  is  the  spectral  absorption  coefficient;  p  =  Cos0;  (b  =  a/(K  +  a);  a  is 
the  spectral  coefficient  of  scattering;  y  is  the  scattering  indicatrix;  Jb  (x)  is  the  spectral  intensity  of  black 
body  radiation.  Intensities  of  radiation  on  boundaries  /^(p)  =  /“(x  =  X//,p),  (p)  =  Jy  (x  =  0,p)  are 

given  for  concrete  calculation  case. 


Figure  3.1 :  Schematic  of  radiation  heat  transfer  in  1 D  geometry 

It  is  assumed  that  radiation  intensity  may  be  presented  in  the  following  form: 

7(X’I (4) 

m= 0  4n 

where  /^(p)  are  the  Legendre  functions. 

Let  scattering  indicatrix  also  may  be  presented  in  the  form 


y (i*»=I(2m+1)yB4(ii)4(4  Yo=1- 

m= 0 


The  Legendre  functions  have  the  ortogonality  property: 


J  Pn(lx)Pm(lx)d^  = 


0, 

2 


2m  +  1 


in  =  n 


(5) 


(6) 
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and  these  functions  are  associated  also  with  the  following  recurrent  formula: 

p  (  l(n)  +  (^  +  1)4+l(n) 

H  mW  2m +  1 

The  Legendre  polynomials  of  low  orders  are: 

7o(p)  =  l,  ^1  (m-)  = 

P4  (p)  =  I^35p4  -  30p2  +3),  P5  (n)  =  i(63p5  - 70p3  +  15(a). 
Substituting  (4)  and  (5)  into  (1),  and  integrating  with  respect  to  total  solid  angle,  one  can  obtain 


V  2m  +  ]  n  (  \ 

E  A  Pm\p) 

m= 0  4k 


dtp m(x)  ,  . 

tl_ar+<p“(,) 


=  i  I  -<■])./,,(  t  )  +  01  y  y  m  pm  ( M  )<P,„  ( c ) 

m= 0  471 


and  taking  into  account  (7) 


(7) 


(8) 


Pm(p)  =  0,  (9) 

where  50  m  is  the  Kronecker  symbol. 

Equation  (9)  must  be  satisfied  at  any  p  ,  therefore 


E 

m= 0 


(m  + 1)  d(Pm+1  (T)  +  m  d9m  1  (2m  + 1)(1  -  ft)ym  )cpm  (t)  -  4ti(I  -  <n)  Jb  (x)50w 


dx 


dx 


(m  + 1) d(Pm^  h)  +  m  d(p^  (T)  +  (2m  + 1)(1  -  &YW  )(Pm  (T)  =  4tt(1  -  m)  Jb  (x)80|W>  m  =  0,l,2,...  (10) 

Unfortunately  system  of  differential  Equations  (10)  is  not  closed.  At  given  m,  each  equation  contains  three 
functions:  cpm_i  (x),  cpm  (x),  cpm+1  (x) . 

For  closing  the  infinite  system  of  equations  it  is  agreed  that  at  m  >  N 

^m±L  =  o,  m  =  N,N  + 1,...,  (11) 

dx 

where  N  is  the  order  of  the  approximation.  The  system  of  Equations  (10)  is  called  as  PN  approximation  of 
the  spherical  harmonics  expansion,  or  simply  PN  approximation. 

Most  simple  Pi  approximation  is  extensively  used  in  radiative  heat  transfer  theory.  In  this  case  ( N  =  1 ) 

d(P!^  +  (1  -  ®Yo  )<Po  (T)  =  47c(1  —  &)Jb  (x) ,  (12) 

ax 

+  3(l-®y1)(Pi(x)  =  0,  (13) 
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where  y0  and  y1  are  the  coefficients  in  expansion  of  scattering  indicatrix 


Y(n'>tO  =  Z  (2m  +  l)ympm  iy)Pm  (f1')  =  1  +  3YllV  • 

m- 0 


(14) 


This  is  so  called  linearly-non-isotropic  indicatrix.  For  isotropic  scattering  ( y0  =  1;  ym  =  0  at  m  >  1 ) 
Equations  (12)  and  (13)  have  the  following  form: 


dcpi  (x) 


dx 


+  (1  -  ro)(p0  (x)  =  4rc(l-  &)Jb  (x) , 


d<Po(Q 


dx 


+  3(Pi(x)-0. 


(15) 

(16) 


If  we  take  into  account  that 


cpo  (x)  =  2n  j  J (x,p)dp  =  G(x)  =  cU (x) , 

-l 

l 

cp1(x)  =  27tj  p/(x,p)dp  =  w(x), 


-l 


because  in  the  Pi  approximation 


1  3 

J  h  0 = CO  *o  (0 + -r<h  (0pi  00  ’ 

4tt  4tt 

^o(0=1’  pi(0=^’ 

then  finally  system  of  equations  for  radiation  transfer  in  plane -parallel  isotropically  scattering  media  takes 
the  form 


dW(x) 


dx 


+  (1  - ®Yo )G(x)  =  4n(l - ®) h  (T)  > 


W(x)  =  - 


1  dG(x) 
3(1-037!)  dx 


(17) 

(18) 


or 


or 


1  d2G(x) 
3(1-007!)  dx2 


+  (1-©Yo)g(0  =  47T(1_®)4(x) 


where  Ub  (x)  =  —Jb. 

c 


Taking  into  account  boundary  conditions  (2)  and  (3),  one  can  get  the  following  boundary-value  problem: 
d2C/(x) 


dx2 


■  =  3(1  -  ®70  )(1  -  a>7i  )U  (x)  -  3(1  -  ra)(l  -  k>7!  )Ub  (x) , 


(19) 
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n  c  dU 
3k  dx 


X=0 


x  =  H:  - 


c  d  U 
3k  dx 


=  2nJ+H --U(x  =  0), 


=  -2kJq  +-U(x  =  H). 


x=H 

where  dx  =  pdx  =  (k  +  a)  dx . 

Let  us  consider  some  numerical  solutions  of  the  boundary-value  problem. 

Figures  3.2-3 .4  show  distributions  of  radiation  energy  density  U,  the  v-projection  of  radiation  flux  W,  and 
divergency  of  radiation  flux  QR  in  plane-parallel  isotropically  scattering  layer  with  temperature 
distribution  shown  in  Figure  3.5. 

Note  that  once  a  boundary-value  problem  (19)  has  been  integrated,  the  v-projection  of  radiation  flux  W  is 
calculated  with  use  of  Equation  (18).  Divergence  of  radiation  flux  can  be  calculated  by  the  following  two 
manners: 


A  formal  differentiation  of  Equation  (18)  leads  to  the  formula 

1  d2G(x) 

V  W  = - - - 

3(l-coy1)  dx2  ’ 

From  Equation  (17)  one  can  get 

V-W  =  4tt(1-©)/4t)-(1-<^0)G(t). 


300000 
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-100000 
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0 


10 
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Figure  3.2:  Volume  density  of  radiative  energy  (1,  U  in  J/cm3),  radiation  flux  (2,  W in  W/cm2),  and  divergency  of 

radiation  flux  (3,  QR  in  W/cm3);  the  optical  depth  x  =  1 
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x,  cm 


Figure  3.3:  Volume  density  of  radiative  energy  (1,  U  in  J/cm3),  radiation  flux  (2,  W in  W/cm2),  and  divergency  of 

radiation  flux  (3,  QR  in  W/cm3);  the  optical  depth  x  =  10 

Figure  3.4  shows  distributions  of  the  functions  for  x  =  1 ,  a  =  10  cm-1. 


x,  cm 


Figure  3.4:  Volume  density  of  radiative  energy  (1,  U  \n  J/cm3),  radiation  flux  (2,  Win  W/cm2),  and  divergency  of 
radiation  flux(3,  QR  in  W/cm3);  the  optical  depth  x  =  1,  a  =  10cM-1 


10-14 


RTO-EN-AVT-1 62 


NATO 

i _ 

OTAN 

Radiation  Modeling  in  Shock-Tubes  and  Entry  Flows 


x,  cm 


Figure  3.5:  Temperature  distribution  in  plane-parallel  layer 
3.1.2  Two-dimensional  radiation  heat  transfer 

Radiation  transfer  is  considered  here  for  axially  symmetric  cylindrical  volume  of  radius  R  and  length  H 
(Figure  3.6).  Such  a  schematic  is  in  common  use  in  different  aerospace  applications,  for  example,  at  study 
of  entry  space  vehicles,  rocket  engines,  plasma  generators,  electro-discharge  devices,  etc. 


Figure  3.6:  Schematic  of  two-dimensional  axially  symmetric  problem 

Radiation  transfer  equation  in  two-dimensional  axially  symmetric  geometry  has  the  following  form: 


r- 


dJ (r, +  1-y2  dJ(r,z,\i, y) 


dr  r  d y 

and  corresponding  to  Figure  3.6  boundary  conditions  are: 


+  pd/(r’^,Y)  +  k/  (r,  z,  p,  y)  =  k Jb(r,z),  (20) 


dz 


z  =  0,  p  > 0 :  / (r,0,p,y)  =  Jq  (r,p,y) , 
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Z  =  H  ,  |Li  <  0  :  J (r,//,p,y)  =  Jq  (r,p,y) , 
r  =  R ,  y <0:  = 

where  p  =  cos0,  y  =  coscp;  /q  (r,p,y),  (r,  jlx, y) ,  7^(z,p,y)  are  the  boundary  intensities. 

Let  us  suppose  that  the  radiation  intensity  J  = /(/?,z,p,y)  can  be  expanded  into  the  spherical  harmonic 
series: 


'=1 

1=0 


1  ^ 

T«/,0^(M)+ZiT(^)(a/.  mcosrncp  +  P/  msinrncp) 

m= 1 


(21) 


where  a/  m  and  P/  m  are  the  unknown  functions  of  spatial  coordinates  r  and  z;  /]m(p)  are  the  associated 
Legendre  functions.  These  functions  have  the  following  orthonormality  relations  [5]. 


2  (/  +  m)! 

2/  +  1  U' 


1  (Z  +  m)! 
2/  + 1  (l-m)\' 


(22) 


where  8//'  is  the  Kronecker  symbol. 

The  following  recursion  relations  are  also  can  be  taken  into  account: 


(2/  +  l)ni;m  (fJ-)  =  (l  —  tn  +  l)//+i  (|J.)  +  (l  +  m) P™x  (^i) , 

0<m<Z-l;  (23) 

(2/ + 1)^7 pr1  (n) = /a  (n)  -  pr-x  M . 

0<m<Z-l;  (24) 

(21  +  l)J77Pim+l (\x)  =  (l  +  m)(l  +  m  + 1  )Pzm!  (|j.)-(Z-m)(Z-m  + 1 )/)'",  (n) , 

0<m<Z  — 1.  (25) 

The  moment  procedure,  which  is  well  known  in  kinetic  theory,  is  used  for  derivation  of  the  PN 
approximation  equations.  An  application  of  integral  operator 


271  1 

WJ{K“a>)K  mcosrncp  +  Pz  msinrncp)dcpdp , 
o  -1 

/=0,1,2,...,  m  =  — /,—/  +  l, ...,0, ...,/— 1,Z,  (26) 

on  Equation  (20)  leads  to  infinite  system  of  ordinary  differential  equations,  which  takes  the  following 
form  at  fixation  of  l  =  N: 
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1  +  8 


m,  0 


2(2/ +1) 


8F, 


m+ 1 


Sr 


-  + 


(m  +  l): 


7m+l 


1  +  8 


m,  0 


2(2/ +  1) 


8F£ f  ( 

— U—  +  (m  +  l)  — — 
Sr  V  }  r 


m+ 1 


+ 


+ 


(Z  +  m  -  l)(Z  +  m)(l  -  8m  o ) 


2(2/ +  1) 


Z7m_1 

- — -  777  ~  1 

Sr  V  }  r 


(Z -m  +  l)(Z -m  +  2)(l-Sm0) 


+ 


2(2Z  +  1) 

Z-m  +  iaF^  Z  +  meFf” 
2Z  +  1  +2Z  +  1  & 


p  17^-1 

Sr  r 


+  K/7/"  =  47riC/68;  08m>0 ,  0  <1<N, 


(27) 


1  +  8 


m,  0 


2(2ZV  +  1) 


ar 


ra+l 

N-l 


dr 


+ 


(m  +  l) 


17/17+1 

P7V-1 


+ 


+ 


(ZV  +  m-l)(ZV  +  m)(l-Sm!0) 


2(2ZV  +  1) 


p  77/77—1  77^77-1 

^VzL_(m_l)^kzL 
ar  V  7  r 


+  N  +  m  dF™_{  +Kfm=(h  t  =  N 


where 


1  2tt 


2ZV  +  1  dz 

Ft m  =  ||  JP/"  (p)cosm(pd(pdp . 

-1  0 

Some  special  cases  are  considered  later. 

1)  PI -approximation  of  radiation  transfer  equation  for  the  plane-parallel  layer. 
In  this  case  N  =  1 ,  and  d/dr  =  0 ,  therefore 
•  1  =  0,  m  =  0 


(28) 


(29) 


dF 


dz 

•  l  =  N  =  l, 


—  +  kF0°  =4nKJb ; 

m  =  0 


1Sf»+KF»  =  0; 


3  dz 

•  l  =  N  =  l, 

2  5+0 


3  az 


m  =  l 
+  kF/  =0; 


In  the  case  under  consideration  must  be  satisfied  conditions  of  axially  symmetry: 

F]  =  Pi  =0, 


because  they  contain  integral  over  even  on  cp  function 
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In  addition 
therefore 

Finally  one  can  get  the  desired  system  of  equation: 

dFr° 


J  /cos(pckp  =  0. 
o 

/)m(|a)  =  0  at  l<m. 


K  =  F()  =  o . 


dz 


+  k/70u  =  AkkJ 


b’ 


kf; 


l^o° 
3  dz 


1  2  71 


where 


F0°  =  j  j  yP0°  (n)d|adcp  =  2n  j  7d(a  =  cU , 


-1  0 
1  2n 


-1 

1 


Fj°  =  |  j  JP°  (n)d|adcp  =  2n  j  y^idju  =  W  . 

-1  0  -1 

Note  that  these  equations  agree  with  Equations  (12)  and  (13)  at  co . 

2)  Pi -approximation  of  radiation  transfer  equation  for  ID  cylindrical  case. 
In  this  case  TV  =  1  ,  and  d/dz  =  0 ,  therefore 

•  /  =  0,  m  =  0 

dFi  F' 


(30) 

(31) 

(32) 

(33) 


+  -L  +  kF0  =  4mc Jb\ 
or  r 


•  l  =  N  =  1, 


m  =  0 


2 

243 


y  dr  r  j 


+  kF{°  =  0 ; 


•  l  =  N  =  1, 


m  =  1 


1 

243 


az?2 


■  2  4 


^5l+25l 

Sr  r 


V 


+  - 


2  dF( 


o 


2-3  dr 


+  KF,1  =  0 


Conditions  of  axial  symmetry  allow  to  get  final  system  of  the  Pi -approximation  for  the  ID  cylindrical 
case: 

^L  +  ^L  +  kF00  =47Oc/,, 
or  r 
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i«£+ri?_0. 

3  dr 


Let  us  take  into  account  that 


1  2  71 


F0°  =  J  J  yP0°  (n)d(id(p  =  2n  J  Jd\x  =  cU , 


therefore 


-1  0 

1  2tt 

f'  = 

JJ 

-1  0 

1 

Lr 

r 

dr 

-l 


(kW)  +  ckU  =  CKUh. 

c  du+w= 0. 


3k  dr 

A  boundary- value  problem  for  ID  cylindrical  volume  with  black  body  surfaces  at  temperature  Tw  is 
formulated  in  the  following  form: 


r  —  0 : 


ld_ 
r  dr 

dU 

dr 


r  dU 


v3k  dr  j 

=  0; 


2)  Prapproximation  for  two-dimensional  cylindrical  geometry. 

In  this  case  TV  =  1 ,  and  with  consideration  of  conditions  of  axially  symmetry,  one  can  derive 

1  drF 1  dF,°  «  „ 

1  -  H - —  +  kF0°  =  4nKJh ; 


(34) 

(35) 

(36) 


r  dr  dz 


o 


Fo=_J_dF}L 
1  3k  dz 


Ft  =  — 


1  S/v 


o 


3k  dr 


Because 


1  2  71  1 

F0°  =  J  J  JPo  (p)dpdcp  =  2tt  J  /dp  =  cU , 
-1  o  -1 

1  271 

F?  =  J  J  JP®  ( (a)  cos  tpdjadcp  =  Wz , 

-1  0 
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1  2tt 


F\  -  |  |  JP'  (p)coscpdpd(p  =  Wr , 


-1  o 


then  the  system  of  equations  takes  more  opaque  for: 


1  d  ,  s  dW 

- ( rWr)-\ - l  +  ckU  =4nKJb, 

r  dr  dz 

3k  dz 


W '=- 


c  dU 
3k  dr 


(37) 

(38) 

(39) 


For  closed  cylindrical  volume  with  black  surfaces  at  temperature  Tw  there  is  the  following  boundary-value 
problem: 


\_d_ 
r  dr 


r  dU 


3k  dr  J  dz 


+  - 


r  dU 
3k  dz 


r  =  0 : 

CD  1  ^ 

II 

o 

t4 

II 

o 

r  =  R: 

z  =  H 

■  k(U  -Ub) . 

2  dU 
3k  dz 

2  dU 
3k  dz 


=  Uh-U. 


(40) 


3)  P3-approximation  for  2D  cylindrical  geometry 
In  this  case  N  =  3  ,  and 


•  1=0,  m  =  0 


•  1  =  1,  m  =  0 


•  1  =  1,  m  =  1 


1  drF{  dF?  0 
- T —  +  — ~  +  kF0  -4nKJb, 

r  or  oz 


1  1  drFj  2  dFo  1  dF?  „  n 

- 2-  + - ^_  + - ^-  +  k F?  =0, 

3  r  dr  3  dz  3  dz 


ac2 


2  3 


dFo  „  F 
dr  r 


1  dF?  1  dF?  1  dF 


0 


.1 


+  — 


3  dr  3  dr  3  dz 


+ - —  +  kF/  =  0 


•  1  =  2 ,  m  =  0 


1  1  drF?  1  1  drF]x  3  dF: f  2  dF |° 


+ ^  + — +  kF2  =0, 


5  r  dr  5  r  dr  5  dz  5  dz 


•  1  =  2,  m  =  1 
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To 


f  dFi 

— — +  2^- 
dr  r 


3  dF,°  3  dF?  2  dF?  3  dF,'  ,  n 

+ - - - —  + - —  + - —  +  kF9  =  0 , 

5  dr  5  dr  5  dz  5  dz 


•  l  —  2  ,  m  —  2 


10 


v  dr  r 


12  dF,1  1  dF}  1  dF?  ,  n 

+ - i - A  + - —  +  kF72  —  0 , 

10  dr  5  dr  5  dz 


•  l  =  N  =  3 ,  m  =  0 


1  1  3  dF° 


1  r  dr  1  dz 


+ - — +  k/v  =  0, 


•  l  =  N  =  3,  m  =  1 


14 


ar2 


2  'A 


dF2  o  F2 
— 2_  +  2^- 

dr  r 


+  - 


12  dF,°  4  dF 


P2-+-^-  +  Kf3‘  =0, 


14  5r  7  & 


•  l  =  N  =  3,  m  =  2 


20 

14 


^  <3r  r 


25t[ff=o, 

1  dz 


•  l =  N  =  3,  m-3 


30 

14 


fSFl  2Fp 

dr  r 


+  kFJ  =0. 


Further  simplification  of  the  system  is  connected  with  the  use  of  properties  of  axial  symmetry. 


3.1.3  Formulation  of  boundary  conditions  for  PN-approximation 

Two  kinds  of  boundary  conditions  are  used  as  a  rule  for  solving  PN-equations. 


The  first  one  is  the  Mark  boundary  conditions  [6].  Mark  suggested  that  boundary  conditions  for  PN- 
approximation  must  be  satisfied  at  angles  corresponded  to  roots  of  the  Legendre  functions: 


WiO=°-  (41) 

As  an  example,  for  the  following  boundary  conditions: 

x  =  0,  p>0:  y(o,p)  =  y+(p), 

t  =  th,  p<0:  y  (xH,p)  =  y“  (p) 

the  Mark  boundary  conditions  are  formulated  in  the  following  form: 

x=o,  p>0:  y(o,p,.)  =  y+(p;), 
t=th,  p<0:  y(Tff,-pi)=y(-pi), 
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where  ju/  are  the  positive  roots  of  Equation  (41). 

For  the  Prapproximation  boundary  conditions  (42)  are  given  at: 

P2(*i)  =  3h?-1  =  0, 

that  is  at  (a,  =  l/sf3 . 

If  the  radiation  transfer  problem  is  solved  relative  to  volume  density  U ,  then 

d  U(t)~ 


t/(x)-|a- 


dx 


and  the  boundary  conditions  are 


x  =  °: 


U(x)~ 


U(  t)+ 


I  dt/(x)~ 
yf3  dx 

1  dU  (x) 


73  dx 

The  second  kind  of  boundary  conditions  was  suggested  by  R.  Marshak  [7]. 
If  physical  boundary  conditions  are  given  in  the  form 

x  =  0 ,  p>0:  /(0,p)  =  /+(p), 

t  =  th,  |i<0:  /(x//,|a)  =  7_(|a), 

then  averaged  moments  of  boundary  intensities  are  calculated  as  following: 

l 

M*  =J/+(|x)[a2Md|x,  |a>0, 


=  J  /  (|a)|a2'  Mia,  |a<0, 


-1 


i  =  l,2,...,(W  +  l)/2. 

For  example,  in  the  Pi -approximation  for  plane-parallel  layer 

dU  (x) 


J  (x,|a)  = 


4n 


[/(x)-|a- 


dx 


and  the  Marshak  boundary  conditions  are  formulated  in  the  form 


M>+  (x  =  0)  =  — 
v  ’An 


1  dt/(t) 


2  dx 


T— 0 


(43) 

(44) 


(45) 

(46) 


(47) 
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M,  (t  =  x„) 


,  1  dC/(x) 

2  dx 


3.1.4  The  Pi-approximation  for  arbitrary  geometry 

In  the  general  case  the  radiation  intensity  is  expanded  into  spherical  harmonics  series: 


(48) 


7v(r,Q)  =  X(2"  +  l)Z 


n- 0 


1 

27t(1  +  50  ^ 

)  (n+  m|)i 

(49) 


where  Y™  are  the  spherical  functions,  which  are  expressed  in  terms  of  associated  Legendre  polynomials 


Yn{^)  =  p<r\v)s™m\ (/,  m>  0, 

y;“m(p,v|/)  =  jp!'")(p)cosmv|/,  m>0, 

— 1  <  p  <  1,  0<\|/<2tt, 

where  80  m  is  the  Kronecker  symbol. 


(50) 


The  Pi-approximation  is  in  common  use  in  the  radiation  gasdynamic.  To  formulate  equations  of  the 
approximation  for  arbitrary  geometry  we  consider  radiation  heat  transfer  equation  for  non-scattering 
medium  in  local  thermodynamic  equilibrium 

Q  •  V/v  (r,f2)  +  kv  (r)  /v  (r,Q)  =  kv  (r)  Jh  v  (r) .  (51) 

Leaving  in  expansion  (59)  only  two  first  terms  one  can  get  the  following  representation  of  intensity 

7V  (r,Q)  =  [cplv  (r)  +  <p2  v  (r)  •  Q]/4t x  ,  (52) 

where  cpl  v(r),  92,v(r)  are  the  unknown  functions. 


Physical  meaning  of  the  functions  is  clearly  defined  at  integration  of  radiation  intensity  over  all  solid  angle: 

J  Jv  (r,Q)df2  =  cpl  v  (r)  =  cUv  (r) ,  (53) 

471 

j  Q Jv  (r,Q)dQ  =  cp2>v  (r)  =  Wv  (r) .  (54) 

471 

So,  in  the  Pi -approximation  radiation  intensity  can  be  presented  in  the  form 

Jv  (r,Q)  =  [< cUv  (r)  +  3WV  (r)n]/47r .  (55) 

The  first  equation  of  the  required  system  is  obtained  by  integration  of  radiation  heat  transfer  Equation  (51) 
over  solid  angle: 

V  •  Wv  (r)  +  kv  (r)cC/v  (r)  =  kv  (r )cUb  v  (r) ,  (56) 
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where  cUh  v  (r)  =  J  Jh  v  (r)dQ . 

471 

The  second  equation  can  be  obtained  after  integration  of  Equation  (51)  with  Q  : 

w”(r)=-5^)vc,»(r^  <57) 

Relation  (57)  shows  that  radiation  flux  is  proportional  to  gradient  of  volume  density  of  radiative  energy. 
Coefficient  D  =  l/3Kv(r)  is  called  as  radiation  diffusion  coefficient. 

Equations  (56)  and  (57)  can  be  presented  in  the  following  form: 

-V  .  VC/v  (r)  +  Kv  (r)ct/v  (r)  =  kv  (r )cUb  (r) 

3k»M  .  (58) 

Typical  boundary  conditions  in  arbitrary  geometry  for  the  Pi -approximation  are: 

1)  No  external  radiation  at  surfaces: 

Jv  (r,Q)  =  0 ,  at  (Q  •  n)  <  0 ,  where  n  is  the  unit  normal  to  surface  s. 

In  this  case 


or 


(r)V£/v(r) 


.S' 


~UA  r). 


2)  Axial  symmetry  for  cylindrical  geometry  or  central  symmetry  for  spherical  geometry: 


(59) 

(60) 


r°Wv(r)  =  0, 

where  r°  is  the  unit  vector  in  radial  direction.  Relation  (61)  can  be  rewritten: 


(61) 


0 

dr 

It  is  of  great  practical  interest  to  consider  equations  of  the  Pi -approximation  in  the  cases  of  the  extreme 
small  and  extreme  large  optical  depths. 

7.  Approximation  of  optically  thin  medium. 

From  exact  equation 


V  •  Wv  (r)  =  Kv  (r  )c[ubv  (r)  -  Uv  (r)] . 
at  xv  □  1  one  can  get  (no  external  radiation) 

VWv(r)  =  cKv(r)C4,v  (r). 


(62) 


(63) 
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Upon  integrating  of  Equation  (63)  with  respect  to  radiation  frequency  we  obtain  the  association  of  integral 
divergency  with  the  Planck  integral  coefficient: 

00 

J  VWv  (r)dv  =  VW(r)  =  4ttkp  (r  )oT4  (r) ,  (64) 

0 


where 


Kt 


■(*)= 


;Jkv  (r  )Ubv  (r)d\ 

0 

00 

1 1 4,v(r)d^dv 


0  4tt 


2.  Approximation  of  optically  thick  medium. 
In  this  case 


00 

\  Kv(r)4,v  (r)dv 
0 

00 

|4v(r)dv 

0 


t/v(r)*t/M(r).  (65) 

Note  that  in  this  case  calculation  of  V  •  Wv  (r)  by  Equation  (56)  may  be  in  error  because  values  Ubv  (r) 
and  Uv  (r )  are  very  close.  In  this  case  one  must  use  the  following  relation: 

w, V(r)=-^v^(r)'  (66) 

Upon  integrating  of  Equation  (76)  with  respect  to  radiation  frequency  we  obtain  the  association  of  integral 
flux  of  radiation  with  the  Rosseland  mean  coefficient 

00  1  C\ 

W(r)  =  (r)dv  = (r)K;‘  (r)' VT(r) ,  (67) 


where 


1  [VKv  (r)][dt/fc;V  (r)/dr]dv 

kr‘  (r) =  ~  x - •  (68) 

\[<WbAT)/dT]dv 

0 

It  should  be  taken  into  account  that  the  Prapproximation  of  spherical  harmonic  method  has  a  poor 
precision  at  small  optical  thicknesses.  To  get  over  the  difficulties,  one  proceeds  as  follows:  the  Pi  method 
is  used  in  that  spectral  regions  where  the  optical  thickness  is  thick,  otherwise  the  emission  approximation 
(radiation  heat  transfer  with  no  absorption)  is  used. 

Examples  of  the  use  of  the  Pi -approximation  for  radiation  heat  transfer  problems  are  presented  in  Part  6.1. 

3.2  The  quadro-moment  method 

The  quadro-moment  method,  suggested  in  [1],  presents  further  development  of  Pi -approximation  for  two- 
dimensional  geometry. 
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Spectral  radiation  intensity  J  =  /(x,r,v)  obeys  the  following  equation  for  non-scattering  medium  under 
LTE  (frequency  dependence  is  left  out  for  charity) 


Y 


v 


8J 
—  + 
dr 


1-y2  dJ^ 
r  dcpy 


SJ  T  T 

+  p - +  KJ  =  KJh  , 

dx 


(69) 


where  p  =  cos0,  y  =  coscp,  Q  =  (cox,<jor),  odx  =  ju,  cor  =  yy  1  - p2  are  the  cosine  of  vectors  Q  with 
respect  to  radial  r  and  axial  x  coordinates;  Jb  is  the  Plank  intensity.  Let  the  problem  is  mirror-symmetrical 
about  the  plane  containing  vectors  x°  and  r°.  It  allows  to  consider  region  for  the  azimuthal  coordinate 
0  <  cp  <  71. 


The  essence  of  the  method  is  that  four  angular  domains  are  introduced  (Figure  3.7)  and  Equation  (69)  is 
integrated  over  these  domains  with  various  weight  functions.  This  is  usually  performed  using  the  moment- 
based  methods  or  multiple  spherical  harmonics  approximations  [3,  8]. 


Figure  3.7:  Schematic  of  the  two-dimensional  axial  symmetry  problem 

In  each  of  these  angular  zones  volumetric  spectral  radiation  densities  are  introduced  by  definition 

U±=Q±{j}  =  2\-^=\j&\i. 


o  4 


I  -  y2  o 


(70) 


The  projections  of  spectral  flux  density  in  radial  direction  are  expressed  as 


W*  =  +R±  {/}  =  2j -p£=J  J^l l-n2dn . 


(71) 


oyi-y  o 


In  the  axial  direction,  the  flux  density  is 


W/  =  Z*  {/}  =  2 j -4M /fid|a  , 


o  \/l- 


■y2  o 


(72) 
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W,±  = 


1  dy  1 


=  -X*  {/}  =  2j  -  l  J /(id|a . 


l-y2  0 


(73) 


Formulas  (70)-(73)  include  indexes  {(~=n, u),  +,  fixing  the  angular  zones: 

(n)-n>0,  (u)-ja<0,  (+)-y>0,  (-)-y<0. 

Consecutive  application  of  operators  1°  Equation  (70)  results  in  the  following  system  of 

equations: 


m  d 
r  dr 


i  aq  /  - 

+  H - h  k 


m  d  (  ~ +  \  m  ~  +  om3 

- M2+n—^  + 

f  of  v  >  r  r)r 


<?M( 


Wr±--Jb 

r  2  ft 


=  0. 


Sr 


w*  -f  4 

z  y 


=  0, 


m  3 
r  3r 


(rW-*) 


SWf  2m 1 


3x 


j/(y  =  0)^^l-n2d|i  +  ^±  )  =  0, 


(74) 

(75) 

(76) 


where 


4±=2j-^=|y(l-^2)dM  ,  M1±=2j1^=jy^^Vd^; 


1  ydy  1 


■\/E 


y2  0 


y2  0 


1  11 —  1  1  1 - 

%  =  2  J  yjl  -  y2  dy  J  /  (l  -  p2 )  dp ,  Mi  =  2  J  /  (y  =  0)  p^/l  -  p2  dp  ; 


=  2j  ydy  J  J\ijl-n2d\i  ,  My  =  2j  dy  J  /  n2d(j. ; 


0  a/i — y2  0 


0  ^i-y2  0 


y>0,  m  =  1 ;  y<0,  m  =  -l;  p>0,  n  —  1;  p<0,  n  —  —\. 


The  resulting  system  of  equations  may  be  used  to  build  a  group  of  calculation  methods.  In  spite  of  the  fact 
that  no  simplifying  assumptions  were  made  for  this  system,  it  can  not  be  considered  as  exact  analogue  of 
the  transfer  equation,  because  the  resulting  system  contains  of  infinite  number  of  ordinary  differentials 
equations. 

Closing  of  the  system  means  that  a  finite  number  of  components  is  used  in  decomposition  of  intensity  into 
series  by  orthonormal  functions.  For  the  full  solid  angle  4tt  the  convenient  orthonormal  basis  is  the  system 
of  spherical  harmonics.  The  first  two  components  of  these  series  (which  corresponds  to  the  Pi  Spherical 
Harmonics  approximation  method)  with  regard  to  evenness  at  cp,  provide: 


/  (  v,  r,  cp,  0)  =  c  (  v,  r  )  +  a  (  x,  r  )  p  +  b  (x,  r 


(77) 
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Approximation  of  radiation  intensity  by  orthogonal  polynomial  at  halved  segments  is  given  as 

/(v,r,cp,0)  =  c(v,r)  +  a(v,r)p|^2p2  -  lj  +  2Z?(v,r)p2^l-p2  ^2y2  -lj  .  (78) 

If  decomposition  similar  to  the  Eq.  (77)  is  used  for  calculation  of  functions  U±,  W±9  Wf,  Mf ,  then 
instead  of  Eqs.  (74)-(76)  one  can  obtain  a  system  of  equation  in  relation  to  coefficients  a ,  b ,  c.  It  is 
reasonably  to  call  this  method  as  a  method  of  quadro-moments,  since  it  is  analogous  to  moment  methods, 
but  corresponds  to  one  quarter  of  a  solid  angle.  Using  the  decomposition  given  by  Eq.  (77),  coefficients  of 
Eq.  (74)  to  (76)  can  be  calculated. 

Let  us  consider  equations  for  the  zero-order  equations  of  the  quadro-moment  method.  Assume 

J±(x,r,cp,0)  =  c±(x,r)  .  (79) 

Using  (79),  Equation  (79)  becomes 


mdrc+  dc+  m  _+  + 

- h  n - c  +2  kc  =  2kJ  h  , 

r  dr  dx  r 

(80) 

T~T±  ~+  xf/i  I  rv+  71 

U  =nc  ,  VU  =  VK..  =  —  c  . 

r  x  2 

(81) 

The  boundary  conditions  are  formed  for  each  of  the  four  directions.  Let  isotropic  radiation  with  intensity 
J°  falls  onto  the  external  surface  of  a  cylinder,  then 

at  r  =  Rc ,  c~  -c~  -  J° ; 
at  x  =  0,  c~  =c+  =J° ; 
at  x  =  Hc,  c~  =c+  =J° . 

At  the  symmetry  axle 


c+  (r  =  ())  =  c  (r  =  0),  c+  (r  =  ())  =  c  (r  =  0)  . 


The  zero-approximation  of  the  quadro-moment  method  (80)  in  relation  to  spectral  functions  has  the 
following  form: 


m  drc+ 
r  dr 


dx 


-~c±  +2kvc±  =  2kvJ 


b,v  ’ 


(82) 


T  T—  ~±  TT/i  it  zi  'tC 

Uv  —  7ZCV  ,  Wv  r  —  Wv  x  —  —  Cv  . 


(83) 


In  group  approximation  (with  the  number  of  groups  Nk  )  index  v  should  be  replaced  by  an  integer  number 
index  k  =  1,2,..., Nk,  which  denotes  a  number  of  spectral  energy  group.  Instead  of  four  Equations  (80)  the 
(4  Nk  )  equations  are  formed.  Integrating  (82)  over  the  whole  spectral  interval  Av  yield 


m  drc+ 
r  dr 


+  n- 


dd± 

dx 


m  „+ 

■ — c~  + 


r 


2 1  kvc+ dv  =  2 1  kvJb  v dv  , 

Av  Av 


(84) 
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where  r  *  -  j"  rv  dv  .  Effective  integral  absorption  coefficient  is  given  as 


Av 


P 


~  + 
c 


J  kvc+ dv . 

Av 


Then  Equation  (84)  becomes 


in  drc+ 
r  dr 


dc+ 

+  n - 

dx 


m 

r 


P  +  2Pc±  =-kpaT4 , 

71 


j  Vfo,vdv 

where  kn  =  ^ -  is  the  average  Planck  coefficient. 

"  \  4,vdv 

Av 


(85) 


(86) 


The  boundary  conditions  for  Equation  (84)  are  obtained  by  integration  of  corresponding  spectral 
conditions  with  respect  to  the  wave  number.  The  advantage  of  such  formulation  of  the  quadro-moment 
method  is  that  unlike  methods  of  averaging  by  the  full  solid  angle  (for  example,  Pi -approximation),  where 
the  average  integral  absorption  coefficient  may  have  gaps  of  the  second  type,  the  given  case  presents  a 
smooth  continuous  symbol  function. 

Examples  of  the  use  of  the  quadro-moment  method  for  radiation  heat  transfer  problems  are  presented  in 
Section  6.2. 

3.3  The  ray-tracing  method 

The  ray-tracing  method  (RTM)  is  usually  used  for  calculation  of  radiation  heat  flux  to  surfaces  of  different 
energetic  devices  and  space  vehicles.  Schematic  of  typical  problem  is  shown  in  Figure  3.8.  In  the  given 
case  the  researched  model  of  Mars  Sample  Return  Orbiter  (MSRO)  has  the  shape  of  a  circular  truncated 
cone  ( a ),  integrated  with  a  spherical  segment  ( b )  on  the  part  of  an  incident  flow,  and  with  the  circular 
cylinder  (c)  in  a  base  area.  For  the  condition  of  axial  symmetry  to  calculate  a  radiation  heat  flux  on  the 
MSRO  surface  it  is  enough  to  perform  calculation  along  a  contour  on  the  surface  shown  in  Figure  3.8. 

Algorithm  of  the  ray-tracing  method  consists  in  the  following.  To  calculate  of  a  radiation  heat  flux  to  any 
element  of  the  MSRO  surface,  the  local  spherical  coordinate  system  with  a  normal  n  is  entered. 
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Figure  3.8:  Schematic  of  research  model  of  Mars  Sampler  Research  Orbiter  Each  ray  Q  is  defined  in  this 
local  coordinate  system  by  the  following  two  angular  coordinates:  the  latitude  angle  0  e  [0,  tt/2]  ,  and  the 
azimuth  angle  cp  e  [0, 2n\ . 

The  radiation  heat  flux  is  defined  by  the  following  formula: 

2tt  7i/2  2tt  1 

J  dcp  J  Jv  (Ry,fi)cos0sin0d0  =  J  dcpj/v  ^RJ9Q^pdp,  (87) 

0  0  0  0 

where  p  =  cos0 ,  Rj  =  ( rj,zj )  is  the  radius-vector  of  the  j-th  point  on  the  MSRO  surface  in  the  laboratory 
co-ordinates,  rj,Zj  are  the  radial  and  axial  coordinates  of  the  j-th  point,  Jv(Rj,Ci)  is  the  spectral  radiation 
intensity  in  the  point  Rj  corresponding  to  direction  Q  . 

Introduction  of  finite-difference  mesh  of  angular  directions  allows  to  integrate  the  spectral  intensity  in 
Equation  (87)  on  angular  variables  and  to  calculate  the  spectral  radiation  flux: 

%-i  Nq-1 

(*/  )  =  Z  (<Pm+l  - <Pm  )  Z  yv  (Qm,«  )(n„  "  1 )  •  (88) 

m= 1  n= 1 

The  direction  cosines  of  the  vector  =(oor)  i  +(oov )  /  +  (oo7)  k  are  calculated  by  the 

m,n  \  xjmn  \  y)mn  V  Zsm,n  J 

following  formulas: 

Mm,n=^-^COSl Pm*  K  )m  n  =  ^/wJsin<P,„,  (©z  ^  =  *1 „  ,  (89) 

where 

Pn  =0.5(mw  +h„+i)  ,  cpm  =0.5(cpm  +cpm+i)  • 

The  z  axis  of  the  local  coordinate  system  of  angular  directions  coincides  with  the  local  normal  to  the 
MSRO  surface. 

To  calculate  spectral  radiation  intensity  a  radiation  heat  transfer  equation  should  be  integrated  through 
inhomogeneous  optical  path.  The  following  solution  of  the  radiation  heat  transfer  equation  can  be  used  for 
this  purpose: 
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s=L 

5 

pm,n)=  f  4v(5)Kv(5)exP 

73 

> 

* 

1 

5—0 

0 

(90) 


where  s  =  0  and  s  =  L  are  the  initial  (on  the  MSRO  surface)  and  final  (on  the  outer  surface  of  the 
calculation  region  or  on  the  MSRO  conjugate  surface)  coordinates  of  the  segment  with  directing  vector 

^ m,n  • 


Finite-difference  scheme  on  s-axis  is  formulated  for  each  ray  Clm  n  .  The  coordinates  of  the  intersections  of 
the  ray  Clm  n  with  all  surfaces  of  the  spatial  finite-difference  mesh  should  be  calculated.  This  is  illustrated 
by  Figure  3.9  (see  points  on  the  ray  a). 


Relations  of  the  analytical  geometry  should  be  used  for  these  purposes.  However,  this  algorithm  is  not 
effective  due  to  very  high  time  consuming,  especially  for  non-structured  meshes.  More  effective  algorithm 
of  the  quasi-random  sampling  was  suggested  in  [1].  A  segment  with  directional  vector  Q  (between  the 
initial  and  final  points)  is  divided  on  Ns  sub  regions  as  it  is  shown  in  Figure  3.9  (see  points  on  ray  b ). 
Then,  the  nearest  node  or  elemental  volume  of  a  spatial  computational  grid  is  searched  for  each  node  sk  . 
In  that  way  temperatures  and  spectral  absorption  coefficients  are  set  to  each  node  of  the  segment 
[.v  =  0,  s  =  L]. 


Figure  3.9:  Two  alternative  finite-difference  schemes  along  s- axis 

Numerical  integration  of  the  solution  (90)  of  radiation  heat  transfer  equation  is  realized  as  follows: 


N,- 1 


Jv  (&m,n  )  “  X  Jb,v,k+ 1/2  {Tk  ~  Tk+ 1 )  ' 


k= 1 


where  J^^^+i/2  the  spectral  intensity  of  the  black  body,  which  is  calculated  using  average  temperature 


of  the  segment  [^,^+1];  Tk  =exp 


sk 

-J  Kv(Y)ds' 


is  the  spectral  transmission  factor. 


Examples  of  the  use  of  the  Pi -approximation  for  radiation  heat  transfer  problems  are  presented  in  Part.  6.1. 
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3.4  Discrete  Ordinates  Method  (DOM) 

The  DOM  on  three-dimensional  rectangular  [9]  and  cylindrical  [10]  computational  grids  was  used 
successfully  for  calculation  radiative  transport  in  a  participating  media.  Ramankutty  and  Crosbie  [11] 
presented  a  modified  discrete  ordinates  solution  of  three-dimensional  radiative  transfer  in  enclosures  with 
localized  boundary  loading  situations.  Another  method,  SN-DOM,  was  used  in  [12]  to  find  the  spectral 
characteristics  of  thermal  radiation  for  a  two-dimensional,  axisymmetric,  free-burning  argon  arc.  Such  a 
method  was  also  used  to  study  radiative  heating  of  internal  surface  of  the  air  and  hydrogen  laser  supported 
plasma  generator  in  [13],  [14]. 

Sakami  et  al.  [15]  proposed  a  modification  of  the  standard  DOM  for  non-orthogonal  computational  grids 
(for  the  tetrahedral  grids).  Further  development  of  the  DOM  algorithm  was  used  to  predict  of  radiative 
heat  fluxes  on  back  surface  of  Mars  Sample  Return  Orbiter  [16]. 

Computational  domain  in  a  rectangular  Cartesian  geometry  is  divided  into  finite  numbers  of  non 
overlapping  tetrahedral  cells.  The  radiation  heat  transfer  equation  in  a  discrete  ordinates  representation  is 
given  by 


r)Tm  r)Tm  r)Tm 

where  I™  is  the  spectral  radiation  intensity,  which  is  a  function  of  position  and  direction  Clm ;  Ibw  is  the 
spectral  blackbody  radiation  at  the  temperature  of  the  medium;  is  the  spectral  absorption  coefficient  of 
the  medium;  r|m,  pm,  are  the  directional  cosines  of  the  direction  Clm ;  w  is  the  wave  number. 

The  condition  at  the  boundary  T  of  computational  domain  is  given  as 


C=  0,  re  r.  (92) 

If  the  surface  of  the  space  vehicle  V  is  assumed  gray  and  emits  and  reflects  diffusely,  then  the  radiative 
boundary  condition  at  the  MSRO  surface  is  given  by 


C=e4r  + 


1-8 


^  HrQm'<0 


X  Om'kfim’lC.  r&T'  . 


(93) 


Spectral  values  are  replaced  by  group  values.  A  control  volume  form  of  Equation  (91)  can  be  obtained  by 
integrating  over  the  tetrahedral  volume  as  follows: 


Z(^«)W=KpVp(/fc>p-/;), 

i= 1 


(94) 


1  f 

where  7™  =  —  I  Imda  is  the  side  area  averaged  group  intensity;  the  number  of  the  tetrahedral  face  is  i 

s,  s. 

and  the  area  of  the  face  is  St ;  7™  =  —  J  7mdo  is  the  cell  volume  averaged  group  intensity;  V p  is  the 

Vp  y 
vp 

volume  of  tetrahedral  cell;  ni  is  the  normal  to  the  side  i  of  the  tetrahedral  cell;  Kp  is  the  group  absorption 
coefficient  in  the  cell;  Ib  p  is  the  group  blackbody  radiation  at  the  temperature  in  the  cell  p  . 
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Solving  Equation  (94)  for  cell  center  group  intensity,  /'”  may  be  evaluated  as 


—  T 

LP  ~Lb,p 


Kpvp  i= t 


(95) 


Therefore,  to  solve  equation  (95),  auxiliary  relations  among  average  intensities  on  sides  of  tetrahedral  cell 
are  necessary.  Tree  cases  arise  with  a  tetrahedral  cell  (Figure  3.10),  depending  on  the  manner  the  cell  is 
approached  by  radiation.  These  cases  were  derived  in  detail  in  [15]. 


If  the  one  face  receives  radiation  from  the  three  other  faces  (Figure  3. 10, a),  then  intensity  on  receiving 
face  4  is  given  as 


For  the  second  case,  when  two  faces  receive  radiation  from  other  two  faces  (Figure  3.10,/?),  intensities  on 
receiving  faces  are 


jin  _ 

h  ~ 


Tm  — 

M  “ 


‘^12  rm  i  “^32  jm 


yS2 


/,  + 


^14  jtn  ,J34  jm 

ySA  1  S.  ' 


2 

\ 
i 

3 


xP+h,p( 

1  %/>)’ 

(97) 

Xp  I b,p  ( 

l~x  p)- 

(98) 

Third  case  describes  the  situation,  when  three  faces  receive  radiation  from  the  fourth  one.  Intensities  on 
receiving  faces  are  given  as 


jin  jin  jin  Tm^  .  T  /i  . .  \ 

h  =I2  =I3  =14%p+Ib,p[l-Xp)- 

Equations  (96)-(99)  include  function  %p  which  can  be  expressed  as 


(99) 
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%P  =■ 


1- 


1  -  e 


p  J 


where  x  is  the  maximum  optical  thickness  in  the  cell  p  indirection 


(100) 


To  predict  the  fluid  flow  and  convective  heat  transfer  in  the  geometrically  complex  systems  of  practical 
interest,  computational  methods  are  normally  implemented  on  non  orthogonal  boundary-fitted  meshes. 
Many  such  problems  also  require  prediction  of  the  radiant  heat  transfer,  where  the  medium  affects  the  heat 
transfer  through  emission  and  absorption.  Radiation  heat  transfer  equation  governs  this  radiant  exchange. 
It  is  desirable  to  solve  this  equation  on  the  same  computational  grid  used  for  the  fluid  flow. 


The  DOM  on  tetrahedral  grids  can  be  applied  to  find  the  radiation  field  in  the  complex  computational 
domain.  In  the  case  of  three  dimensional  geometry  octahedral  cell  is  divided  into  six  tetrahedrons  ABCF , 
AFDE ,  FCDA ,  DEGH  ,  GEFD ,  GFCD  (Figure  3.11,  a).  The  cylindrical  meshes  also  have  hexahedral 
cells  adjoining  z  axis  (Figure  3.12,  a).  Hexahedrons  are  divided  into  three  tetrahedral  cells  ABCD , 
BCEF ,  CDEB  (Figure  3.11,  b).  Using  quasi  uniform  meshes  (Figure  3.12,  b)  is  preferable  in  the  case  of 
axial  symmetry. 


a  b 

Figure  3.11 :  Tetrahedral  cell  partition 


a  b 

Figure  3.12:  Cross  section  of  cylindrical  computational  domain  that  is  perpendicular  to  the  z  axis  (a).  Cross 
section  of  quasi  uniform  grid  constructed  in  cylindrical  computational  domain 


The  accuracy  of  the  discrete  ordinates  solution  depends  on  the  choice  of  the  quadrature  scheme.  Although 
this  choice,  in  principle,  arbitrary,  a  completely  symmetric  quadrature  is  preferred  in  order  to  preserve 
geometric  invariance  of  the  solution.  The  quadrature  schemes  used  in  [16]  are  based  on  the  “moment¬ 
matching”  technique,  whereby  the  ordinates  are  chosen  so  as  to  integrate  as  many  moments  of  intensity 
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distribution  as  possible.  Although  any  type  of  the  quadrature  can  be  applied  the  SN  quadrature  sets  were 
used  in  [16].  The  recognition  algorithm  was  created  to  identify  the  type  of  Clm  orientation  in  the  cell.  This 
algorithm  permits  to  apply  appropriate  characteristic  Equations  (96)-(99)  for  finding  unknown  average 
face  intensity.  If  the  solver  visits  the  cells  in  the  correct  order,  all  terms  on  the  right  side  of  the  Equation 
(95)  will  be  known  so  center  intensity  can  be  found  by  direct  substitution.  This  allows  the  solution  to  be 
obtained  by  moving  from  cell  to  cell  in  the  optimal  “matching  order”. 

For  the  irregular  grid  the  order  is  not  obvious.  To  find  a  solution  on  irregular  grids  the  equation  could  be 
solved  by  repeatedly  sweeping  across  the  grid  without  regard  to  the  optimal  order.  Creating  of  a 
“marching  order  map”  gives  the  optimal  order  in  which  the  cells  will  be  visited.  Once  the  quadrature 
scheme  is  fixed,  a  marching  order  map  can  be  constructed  for  each  intensity  direction,  and  this  map  can  be 
written  to  file  for  repeated  use.  To  create  the  map,  the  boundary  cells  are  swept  first  to  find  the  starting 
location.  The  map  is  easy  to  generate,  it  needs  to  be  revised  only  when  the  spatial  grid  or  quadrature 
scheme  have  been  changed.  Because  it  reduces  solution  costs,  the  creation  of  the  map  is  the  preferred 
alternative. 

Examples  of  the  use  of  the  Pi -approximation  for  radiation  heat  transfer  problems  are  presented  in  Section 
6.4. 

4.0  THE  RANDOM  MODELS 

4.1  Formulation  of  Random  Models  for  Atomic  Lines 

Random  model  of  molecular  lines  was  originally  offered  and  investigated  in  the  Refs.  [17,18].  Detailed 
description  of  the  random  approach  is  presented  in  [19].  Being  based  on  the  specified  references,  we  shall 
obtain  formulas  of  generalized  random  model  of  rotational  molecular  lines  and  spectral  atomic  lines.  Note 
that  the  distinctive  property  of  band  of  atomic  lines  is  that  the  half-widths  of  the  atomic  lines  differ 
significantly.  In  the  case  of  rotational  molecular  lines,  theirs  half- widths  differ  insignificantly. 

Let  a  spectral  range  Aco  contains  N  atomic  lines.  Then  spectral  transmission  along  the  homogeneous 
optical  path  L  is  calculated  as  follows: 

N 

T((o)  =  YlT((o,(oi),  (1) 

i= 1 

where  r(co,co/)  =  exp[-K/(o3,(jo/)L],  k^cd,^)  is  the  spectral  absorption  coefficient  of  /- th  line,  co,  is 
the  wavenumber  of  i- th  line. 


Integrated  transmission  in  the  spectral  range  Aco  is  determined  under  the  formula 


T=  J  r(co)dco 


Aco 


and  average  transmission  is  determined  as  follows 

Tac°=~L)  J  r(w)da)- 


Aco 


(2) 


(3) 
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It  is  obvious,  that  average  transmission  (3)  will  depend  on  how  the  lines  are  placed  inside  the  range  Aco . 
Intuitively  also  it  is  clear  that  at  large  number  N  this  dependence  is  not  strong.  The  truth  of  such  an 
assumption  easily  to  check  up  by  means  of  line-by-line  integration  for  any  set  of  lines.  This  fact  allows 
give  a  probabilistic  interpretation  for  averaged  transmission.  For  this  purpose  it  is  necessary  to  formulate 
probability  of  realization  of  just  given  configuration  of  lines  inside  the  range  Aco ,  and  then  to  estimate 
mathematical  expectation  of  averaged  transmission  for  every  possible  configurations. 

Let  us  define  probability  for  i- th  line  to  have  intensity  in  the  range  dSt  in  the  vicinity  of  St ,  half-width  in 
the  range  dyz  in  the  vicinity  of  yt ,  and  location  in  the  limits  d cot  in  the  vicinity  of  co,  as  follows: 

d Pt  =  dPi((oi,Si,yi)  =  pi((oi)pi(Si)pi{yi)d(oidSidyi,  (4) 

where  statistical  independence  of  realization  of  the  location,  intensity  and  half-width  of  the  line  is 
supposed,  /?(co/),/?(5/),/?(yl-)  are  the  probability  densities  for  realization  of  the  location,  intensity  and 
half-width  of  the  line.  Then  the  mathematical  expectation  of  the  integrated  transmission  can  be  determined 
as  follows: 


T=\ 


Aco 


A®!  A  A,  Ay,  Aco^  ASN  Ay  N 


I  I  I  -  I  I  I  |r(®)nA(®;)d®;ft(^)d^A(Yi)dYi 

i= 1 


>dco . 


(5) 


Now,  let  us  assume  that  the  probability  of  realization  of  the  parameters  of  each  line  ( CO; ,  S[ ,  y, )  do  not 
depend  on  realization  of  the  parameters  of  any  other  line.  It  allows  simplify  last  relation: 


f=j< 

Aco 


N 


XT  Iff  T  (  ’  S‘  ’  Yi )  Pi  ( )  P  ( Si )  P  ( Y, ) d®;  dSt  dY , 

i= 1  Aco,  AA,  Ay, 


dco, 


(6) 


where  the  explicit  dependence  from  St  and  yt  is  entered  into  the  list  of  arguments  of  spectral 
transmission. 


Among  lines  located  in  the  spectral  range  can  be  as  strongly  differing  lines,  and  the  lines  with  rather  close 
parameters.  Let  us  enter  conditional  division  of  the  lines  into  groups,  in  limits  of  which  the  probabilities  of 
realization  of  the  lines  parameters  can  be  described  by  the  uniform  law.  Thus  we  shall  remind,  that  in  the 
limits  of  any  allocated  group  all  lines  are  indiscernible  (the  statistical  independence  of  the  lines  parameters 
already  was  used  at  transition  from  Equation  (5)  to  Equation  (6)).  Then,  Equation  (6)  can  be  rewritten 

(  G  Ng  ] 

nilj  j  j  T((o,co;-,5,;-,Y/)A(®/)A(^)p;(Yi)d®id‘S'idY;-  fd®,  (7) 

Aco  18= 1  i= 1  Aco  i  A  A,  Ay,-  J 

where  Ng  is  the  number  of  lines  in  g- th  group,  G  is  the  number  of  the  groups. 


The  further  analysis  of  Equation  (7)  will  be  connected  with  concrete  definition  of  the  probability  densities. 
We  shall  accept,  that  all  spectral  ranges  Aco t  are  equal  among  themselves  Aco t  =  Aco, i  =  l,2,..., Ng  thus 
they  are  distributed  equiprobably,  that  is  pt  (co, )  =  1/ Aco ,  then 

f=Jn^J  J  J  T(a,a',S',y')p(S')p(y')doi'dS'dy' 

Acog=l  Att>ASgAyg 

where  A Sg  and  Ay g  are  the  ranges  of  change  of  the  intensity  and  half- width  in  limits  of  the  g- th  group. 


g 

dco ,  (8) 
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Expression  in  square  brackets  has  meaning  of  the  average  transmittance  in  spectral  group 

~\N, 


Aco 


kj  j  J  T (©,©', S', y') p(S')/>(y')d(o'dS'dy' 

Aco  A Sg  Ayg 

Let  us  enter  a  new  scale  of  wave  numbers,  connected  to  centre  of  the  spectral  range 

Aco  —  comax  —  CDmjn  .  V  —  CO  —  , 

where  co0  =  0.5 (co^  +  comax ) ,  comin,comax  are  the  boundaries  of  the  spectral  region.  Then 

Vmin  =  -0.5  Aco  ,  vmax  =  +0.5  Aco . 

If  the  average  distance  between  lines  in  g- th  group  to  define  under  the  formula 

8s=4co/Ws, 

then  vmjn  —  — 1 0.5Ng5g ,  vmax  —  +0.5A^Sg  . 

Now  it  is  possible  to  copy  (9)  in  the  following  form: 


7k(co)  = 


V'5 


k  J  J  T.  J  r(o3,o3',5',Y')p(5')p(Y')dv'dS'dY' 


M  J  J  § 

g*sg  Ays°8  JLv8 


j 

1-  |  J  T(a>,ai',S',y')p(S')p(y')dS'&r' 

dv'  > 

A Sg  Ayg 

N„ 


J _ 1_  z 

*+§7  1 


At  rather  large  number  Ng  it  is  possible  to  pass  to  the  exponential  form  in  the  Equation  (13): 

rx(v)=exp(-iys/5s), 

where 


h — N080 
2 


w8(v)=  i 


-5“AL 


1-  J  J  r(co,co',S',Y')/?(S')/?(Y')dS'dY' 


dv'  : 


1-  j  j  r(co,co',5',YX5')^(y')d5'dY' 

AS„  Ay„ 


dv'. 


(9) 


(10) 


(11) 


(12) 


(13) 


(14) 


where  Wg  is  the  equivalent  width  of  lines  of  g-th  group.  The  transition  from  final  limits  of  integration  in 
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(14)  to  infinite  can  be  carried  out  formally,  by  requiring  that  Ng  — »  oo  . 

However,  it  is  enough,  that  the  following  condition  should  be  satisfied 

yO  Aco.  (15) 

Therefore  such  transition  can  be  made  basically  even  at  Ng  =  1 . 

Above  already  was  spoken  that  it  is  possible  to  divide  all  lines  in  the  range  Aco  on  number  of  groups.  For 
example,  it  is  possible  to  use  the  assumption  of  affinity  of  half-width  of  lines  in  limits  of  each  group  and 
about  identical  probability  law  of  lines  intensity  distribution: 

p(Y')  =  §(y'-yg),  (16) 

p{S  )  -  Sgl  exp(-S'/Sg )  •  (17) 

Condition  (16)  means  that  inside  given  group  all  lines  have  identical  half- widths.  The  relation  (17)  is 
frequently  used  in  the  theory  of  random  modeling  of  rotational  lines  of  a  molecular  spectrum  and 
physically  is  substantiated  in  Ref.  [17].  With  regard  to  (16)  and  (17)  one  can  obtain  instead  of  Equation 
(14)  the  following  relation: 


+00 

«+)=  \ 


00  . 

l-Jr(y,y',S',yg)— exp 


f 

v  V 


dS' 


dv'. 


(18) 


0  8 

The  intensity  of  a  line  S,  by  definition,  is  integrated  absorption  coefficient  in  this  line,  therefore 

k(  v,  v',  S',yg^  =  S'f  (v,  v',yg  ) , 

+00 

where  the  function  of  the  line  contour  satisfies  to  the  following  condition  J  /’(v.v )dv'  =  1  . 


(19) 


Therefore,  by  substituting  (19)  in  (18),  one  can  obtain: 


f  S' ^ 


«+)=  J  U 57eXp[  - 57  |D -exj+sr (V.V'.T«  ))+s'  |dv'  =  f  J 


V  8  J 


r  V( 

v’v'’YA 

L 

b+V 

(v’v'’Yg 

)L 

dv' .  (20) 


The  further  transformation  (20)  is  possible  at  detailing  of  spectral  dependence  of  absorption  coefficient  in 
atomic  lines.  Exact  and  approximation  expressions  for  equivalent  width  of  lines  of  Lorentzian,  Doppler, 
and  Voigt  contours  were  obtained  in  the  theory  of  statistical  modelling  of  spectral  lines  [17,18]. 

Thus,  the  mathematical  expectation  of  integrated  transmission  in  a  range  Aco  takes  the  following  form: 

(21) 


T  =Aco]^[exp 

8= 1 


v 
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Hence,  average  transmission  in  any  range  does  not  depend  on  wave  number 


z  Aco 


TIexP 

g=i 


IT'. 

£=1 


(22) 


The  given  formulation  of  random  model  contains  a  large  variety  of  random  models:  from  simple  random 
model  (G  =  1)  up  to  generalized  random  model  (G  =  N)  of  superposition  of  lines  of  various  bands.  In 
conformity  to  the  formulation  of  generalized  random  model,  given  by  Goody  [17],  now  it  is  possible  to 
assume,  that  the  size  of  the  spectral  range,  in  which  the  average  transmission  is  sought,  coincides  with 
average  distance  between  lines  in  group.  It  means  that  the  random  model  may  be  used  even  for  one  line. 
In  other  words,  it  is  possible  to  determine  average  transmission  in  a  spectral  range,  containing  only  one 
line,  being  the  member  of  a  spectral  band  with  average  distance  Aco  between  lines. 

4.2  Numerical  Simulation  Method  for  Calculation  of  Radiative  Heat  Transfer 
in  Plane-Parallel  Non-Uniform  Layers 

Now,  let  us  consider  how  it  is  possible  to  use  the  generalized  random  model  for  calculation  radiative  heat 
transfer  in  plane-parallel  non-uniform  layers  of  low-temperature  plasma.  For  a  basis  let  us  consider  the 
equation  for  radiation  transfer  in  the  following  form: 

=  k<oK» -ya)>  (23) 

where  are  the  spectral  intensity  of  radiation  of  the  plasma  and  absolutely  black  body;  is 

spectral  volumetric  absorption  coefficient;  a  is  a  direction  of  radiation  distribution. 

For  determination  of  the  spectral  and  integrated  radiating  characteristics  in  a  plane-parallel  non-uniform 
layer  the  method  of  half-moment  is  used  [8].  The  Equation  (23)  can  be  approximated  by  equivalent  system 
of  the  differential  equations  for  spectral  half-moment  characteristics  : 


l,co 


=  -M 


n,  co 


2n 

+ - 

n  + 1 


J, 


b,  co’ 


n  =  1,2,..., 


where 


dM 


n+ 1,(0 


dx„ 


=-m: 


,  M-}T  j 

'  n  + 1 


n  =  1,2,..., 


(j,  =  cos0,  dxC0=Ka)dx;  =  2nMn^  =  2ti(-1)"  ; 

o 

_f/+,  0  <  0  <  rc/2, 

^  CO  I 

[Ja,  7r/2  <  0  <  tc; 


0  is  the  angle  between  direction  of  radiation  distribution  and  v-axis. 


Let  us  assume  that  external  radiation  is  characterized  by  the  following  spectral  intensity: 


(24) 

(25) 


(26) 
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x  =  0,  /ro(|^,x  =  0)  =  /+(x  =  0),  |-i  <=[0,1], 

x  =  H,  Ja(ii,x  =  H)  =  J~(x  =  H),  ju  <=  [—1,0], 
where  His  the  thickness  of  the  layer. 

To  limit  number  of  the  solved  equations  in  system  (24)  it  is  necessary  to  use  any  closing  condition,  by 
putting  all  half-moment  functions  equal  to  zero  since  some  number.  Such  closing  condition  is  possible  to 
obtain  approximating  an  angular  dependence  of  intensity  by  a  series  on  cosines  with  limited  number  N  of 
terms.  Having  integrated  this  series  over  two  half-spherical  angles  (0<0<7i/2,  and  n/2<6<n)  with 
different  weight  functions  \in ,  it  is  possible  to  obtain  the  following  closing  relations: 


N  a+,  0(0  ) 


V,Q3  (^m)  ~  \r  ■  i 

P  N-j+n+l 

n  (-l)N-j+na~(x  ) 

l)  , 

P  N-j+n+l 


(27) 


(28) 


where  a~j  (xro)  and  aj  (x^)  are  the  coefficients  of  the  decomposition.  By  excluding  them  from  system 
(27),  (28),  it  is  possible  to  find  connection  between  half-moment  characteristics: 

M2, co  =~~M0xo  +  M2, co  =_5^()’C° 

In  the  second  approximation  of  the  half-moment  method  ( N  =  2,  n  =  0,1,2)  the  systems  of  equations  has 
the  following  solution: 


2  x 


M0,a>  (4  =  £  Pm  f  4,co  (0 1 Kco  ( V)  exP 

m-\  o 


2  x 


Mu„  (4 = Z 1 4,co  b'K  (^')exP 


m= 1 0 


2  H 


M0,co  (4  =  £  Pm  1  4,  co  b')Kco  (x')exp 


m— 1 


2  H 


(X)  =  “Z  J  4,03  (x')K«>  (x')exp 

where 


m= 1 


f  ^ro  (-X")<±V" 

x' 

x' 

-PmfKco(V)dxff 

X 

-Pm]Km(x'')cb:" 


d^'  +  ZPmCexp 

m=l 

2 

dv'  +  ^  C  exp 

m=l 

2 

<**'  -  Z  PmCm  eXP 

m=l 
2 

dr'  +  ^  C“  exp 


“Pm  jK03  (x')dr'  ,  (29) 

0 

Pm}Ka,(x')dx'  ,  (30) 

H 

Pm  J  Kco  (x')dr'  ,  (31) 

X 

H 

-PmjKro(x')dr'  ,  (32) 


4,03  ^4,03  (4  5  Pi 

c+  K  ~  P24+  c+  4+  ~Pi4+ 

1  P2-P1  ’  2  P2-P1 


3-73,  P2  =  3  +  73  , 


*b~  +  P2*r  r-  4T  +  P7T 
1  P2-P1  ’  2  P2-P1 
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Fq  =  2nj /+  (jc  =  0)dp,  F{+  =  2tt|  J+  (x  =  0)pdp , 

o  o 

1  l 

Fq  =  2ttJ  J-  (x  =  //) d|Li,  F{  =  2ttJ (x  =  //)pdp  . 


o  o 

Let  us  now  find  integrated  half-moment  characteristics  in  the  spectral  range  Aco  .  Calculations  for  the  four 
half-moment  functions  (29)-(32)  are  identical,  therefore  we  shall  obtain  a  required  relation  only  for 


KM- 


Integrating  of  Equation  (30)  in  the  limits  of  the  Aco  one  can  derive: 


2  x 


KM=  1  K(M=1S  j '  bAx’)KJx’M 

m= 1 0  Aco 


Aco 


dcodx'  + 


+Z  |  Cexp 

m=\  Aco 


dco. 


The  first  term  in  the  right  hand  side  of  Equation  (33)  correspond  to  own  radiation  of  the  medium: 


A  =f  f  4,co(x')Kco(x')exP 


0  Aco 


-Pmj’Kco(x'')(k" 


1  AWm(x\x) 

dodx  =-|4(0(x)- - A_ — K  , 

Pm  ^ 


where 


W 


( x\x)=  J  \  1-exp 


Aco  | 


-Pm}Kco(*'')d*'' 


^dco, 


IbAx)=-Z-  1  4,co(*)d®- 

Aco  f 

Aco 


(33) 


(34) 


(35) 


At  calculation  of  integral  Ix ,  conveniently  to  use  the  following  formula 


Pm  i—\ 


where  /  is  the  number  of  a  mesh  nodes. 


(36) 


Let  us  present  the  absorption  coefficient  k00(jc)  as  the  sum  of  absorption  coefficient  in  lines  and 

one  in  a  continuous  spectrum  k^(x).  It  is  assumed,  that  (a;)  =  kc  (a:)  ,  that  is  it  does  not  depend  on 
wave  number  in  the  limits  of  the  Aco .  Then  the  expression  for  equivalent  width  (35)  can  be  presented  in 
the  following  form: 
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Wm(xi’xf): 


1-exp 


*/ 

P  -Pm  J  <,(*')dx' 

>  Aco  +  exp 

1 

1 

"CO 

s 

7\ 

-- 

Xi 

II 

1-exp 

xf 

"Pm  j  Ki: 

((x')dx' 

Aco 

Xi 

K(*i>xf)>  (37) 


dco .  (38) 

For  calculation  of  the  equivalent  width  of  lines  W*n  (  ,v,  .  )  one  can  use  the  group  random  model: 

K(xi'Xf)  =  tiWw(xi,xf),  (39) 

S=1 

where  W*m  (xz-,xy )  is  the  equivalent  width  of  the  g-th  group;  G  is  number  of  atomic  lines  groups. 


Now  consider  transformation  of  integral  corresponding  to  the  contribution  of  external  radiation.  If 
intensity  of  external  radiation  poorly  changes  in  the  spectral  region  Aco ,  the  factors  c*  may  be  considered 
as  constants,  then 


j  exp 

i 

1 

"CO 

3 

7\ 

S 

X 

_ I _ i 

Aco 

0 

dco. 


Because  (V)  =  (V)  +  (x')  ,  then 


1 2  =exp 

Xf 

-Pm  \  <,(*0^' 

•  j  exp 

xf 

-Pm  J 

0 

Aco 

0 

dco . 


In  this  formula  integration  over  wavenumber  gives  total  transmission  in  the  region  Aco,  that  is  in 
conformity  with  (21) 


Then 


G 

T^xx  =  0,Xy  )  =  Aco]^[exp 
g= 1 


Wg,m(* 1  =  °.*/) 
Aco 


1 2  =  exp 


"Pm  f 


TIexp 
8= 1 


Wg  m  ^  X^  0,  Xy  ^ 


Aco 


(40) 


Thus,  the  required  formula  for  determination  of  the  half-moment  characteristic  MyA(0(xy)  has  the 
following  form: 


m= 1  Pm  i= 1 


+ 


+Zexp 

m= 1 


“Pm  J 


G 

n«p 

8=  1 


Wg,m(x\  =  °’Xf) 


Aco 


Aco, 


(41) 
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where  the  equivalent  width  may  be  calculated  by  the  following  approximate  formulas  [19]: 

wm (xf,Xj )  =  AAam Aeo  =  [l - TAmm (xf,Xj )] Aco  =  1 1  -Ylexv[~Dk  (xf’xj )] | Aco > 


k= 1 


0fh/.v)=,tof+fef  + 


ryn  jyn 
uL,k1JD,k 

Kjm,k 


-  ^Lk  -  D™k(xf  ,Xj)~- 


Xm,k 


D D,k  -  D D,k  (xf  ’  x  j  )  ~~  '  G aD,k 


In 


1  + 


X2 

%ra,& 


B  ' 1 

»  Xm,k  =1m,k(xf’Xj)  =  ^rL\ 


Atomic  parameters  of  the  model  are  calculated  as  follows: 


m 

aL,k 


[T(x)]GD,k[T(x)]  dx. 


j= i 


Nt 


GD,k{T)  =  ^Dj  =3.58x10 

7=1 


'Mr  ■ 


where  Na ,  Ne  are  the  concentrations  of  atoms  and  electrons;  fj  is  the  oscillator  strength;  gj  is  the 

statistical  weight  of  the  level  with  energy  Ej  ;  C^^C\  j  are  the  Stark  constants  for  atomic  ( a )  and  ionic 
(/)  lines;  co 0j  is  wavenumber  of  the  atomic  line  centre;  M  is  atomic  weigh;  k  is  the  Boltzmann  constant;  Q 
is  the  total  partition  function. 


It  must  be  stressed,  that  the  group  functions  depend  only  on  temperature  and  do  not  depend  from 
population  of  absorbing  energy  levels.  It  means  that  at  calculation  of  radiation  transfer  they  can  be 
determined  at  once  beforehand.  Using  of  the  group  functions  allows  reach  appreciable  economy  of 
computer  resources. 

Examples  of  the  use  of  the  random  models  of  atomic  lines  for  radiation  heat  transfer  problems  are 
presented  in  Section  6.5. 
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4.3  The  Macro-Random  Model  for  Describing  of  Radiative  Heat  Transfer  with 
Vibrational  Band  Structure 

The  general  idea  of  the  Macro-Random  model  is  that  each  vibrational  band  is  considered  as  an  isolated 
line  of  absorption.  The  average  transition  in  a  spectral  region  is  obtained  by  multiplication  of  separate 
band  transitions.  The  macro-random  model  is  computationally  very  effective  as  it  allows  determine 
averaged  radiative  characteristics  without  the  use  of  " line-by-line "  integration  or  any  group  models.  The 
calculations  are  performed  as  if  there  is  a  single  spectral  region.  For  calculation  of  radiative  heat  transfer 
in  non-uniform  volumes  the  Curtis  -  Godson  method  is  used. 

Basic  peculiarity  of  the  Macro-Random  model  is  the  fact  that  the  model  is  formulated  with  reference  to 
one  line  within  a  spectral  range  Aco .  Such  a  random  model  was  introduced  and  justified  above.  Any 
contradiction  with  the  traditional  formulation  of  random  models  is  not  present,  since  in  the  given  random 
model  it  is  supposed,  that  in  a  considered  spectral  range  Aco  only  one  spectral  line  has  got  from  some 
large  population,  average  distance  between  lines  in  which  equally  Aco . 

The  opportunity  of  application  of  such  a  model  to  calculation  of  averaged  spectral  characteristics  was 
checked  by  comparison  with  results  of  “line-by-line”  calculations  [20,21]. 

As  it  was  mentioned  above,  the  basic  idea  of  the  macro-random  model  is  that  the  vibrational  band 
averaged  on  rotational  structure  is  considered  as  an  isolated  line  of  absorption.  To  use  the  generalized 
random  model  it  is  necessary  to  set  a  contour  function  of  the  absorption  band  /(v,v',S,y') ,  where  the 
intensity  S  and  the  half-width  y  are  parameters,  specifying  the  integrated  absorption  coefficient  and  half¬ 
width  of  the  band.  It  is  necessary  also  to  fulfil  the  condition  that  the  half-width  of  the  band  should  be 
much  less,  than  size  of  the  spectral  range  Aco  . 

In  the  study  the  approximation  of  the  function  of  vibrational  bands  absorption  offered  by  Edwards  [22]  is 
used: 


Ci 


k(v,v')  =  —  exp 


Ca 


V-V 


Ca 


(42) 


where  v'  is  the  centre  of  the  band  of  absorption,  C\  and  C2  are  parameters,  determining  intensity  and  half¬ 
width  of  the  band.  Let  us  assume  that  v'  =  0,  then  easily  to  show  that 


00  ^ 

s  =  2j  exp 


oC3 


f  \ 

V 

V  c3  J 


dv  =  2C, . 


(43) 


The  approximation  (42)  allows  to  accept  that  y '  =  C3 ,  then 


(44) 

If  to  assume,  that  the  probability  density  of  intensities  distribution  of  vibrational  bands  has  the  same  kind, 
as  well  as  for  rotational  lines  in  vibrational  bands,  i.e. 


S' 

k(v’v')=— exp 


^(5')=^-exp 


\  SSJ 


then  for  obtaining  the  required  equivalent  width  it  is  necessary  to  calculate  integral  (18),  by  substituting  in 
it  the  following  contour  function 
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Av’v'^)  =  ^ex  P 


v-v 


Then 


^=2ygln 


\  SgL^ 

1+  g 


2y, 


(45) 


'sj 


Thus,  average  transmission  of  the  spectral  range  Aco ,  containing  vibrational  bands  can  be  calculated  under 
the  formula 


-  G  f 
t&o> = n-p 

8= 1  V 


w  ' 

8 

Aco 


where  the  equivalent  width  of  the  g- th  vibrational  band  Wg  can  be  calculated  under  the  formula  (45). 


To  calculate  an  average  transmission  of  a  non-uniform  optical  path  the  Curtis  -  Godson  method  is  used 
[23,24].  The  general  idea  of  this  method  consists  in  introduction  for  each  line  (or  each  band)  averaged 
parameters  on  the  optical  path:  average  intensity  of  the  line  Sg  and  the  average  half-width  yg  . 

The  introduction  of  averaged  parameters  Sg  and  yg  actually  means  transition  to  consideration  of  some 
effective  optical  path.  As  it  was  offered  by  Curtis  and  Godson,  to  find  averaged  parameters  it  is  necessary 
to  equate  the  integrated  absorption  (or,  that  same,  equivalent  width)  of  the  effective  optical  path  and  one 
for  the  true  non-uniform  path  in  two  limiting  cases:  in  the  cases  of  optically  thin  and  optically  thick  layers. 


Let  us  consider  obtaining  of  the  Curtis  -  Godson  relations  with  reference  to  one  isolated  line.  An  exact 
expression  for  the  equivalent  width  has  the  following  form: 


+00 

WT=  J  <j  1-exp 


Idv  =2Nl-exp 


-J  K(v,S(x),y(x))ck 

L 

Then,  for  the  effective  homogeneous  path  one  can  write 

00 

w  =  2  j  {l-exp[-K(v,5,y)L]}dv . 
0 

In  either  case,  the  wave  number  v  is  measured  from  centre  of  the  line. 


J  K(v,5(x),y(x))dv 


>dv  (46) 


(47) 


For  obtaining  of  the  averaged  parameters  S  and  y  the  following  condition  WT  «  W  must  be  satisfied,  or 


J  exp[-ic(v,S,y)L]dv  «  j  exp 
o  o 


-  J  k  (v,  S  (x),  y(x) )  dv 

L 


dv. 


(48) 


It  is  obvious,  that  from  this  relation  it  is  possible  to  determine  required  values  only  approximately.  Let  us 
consider  the  following  two  opposite  approximations:  the  optically  thin  (k(v,5,y)LD  l)  and  the  optically 
thick  (k(v,5,y)lB!  l)  paths. 
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The  equivalent  formulation  for  these  approximations  looks  as  follows: 

jK(v,5,(x),y(x))drD  1  foranyv,and  Jk(v,S(v),y(v))cUD  l,atv<v1. 

L  L 

The  sense  of  the  last  inequality  is  that  one  can  find  the  wave  number  Vj  □  max{y(v)},  at  which  the 
optical  thickness  of  the  volume  will  become  less  unit.  In  this  case  the  radiation  transfer  must  be 
investigated  only  in  far  wings  of  the  line  ( v  >  Vi ). 

Let  us  consider  the  following  integral 

exp  [-k(  v,  5 ,  y )  l]  -  exp 

o  l 

In  the  case  of  optically  thin  path  one  can  obtain 


-|K(v,5'(x),y(x))dx 


>dv  =  0. 


(49) 


Lj  y)dv  =  JdxJ  ic(v,S(x),y(x))dv 

0  L  0 

or  using  the  definition  of  intensity  (19): 

S  =— Js(x)dx.  (50) 

Ll 

In  the  case  of  optically  thick  path  (at  v  <  Vi )  the  total  absorption  of  radiation  will  be  observed,  so  both 
exponents  in  (49)  are  practically  equal  to  zero  and  it  is  possible  to  write  the  following  approximate  relation: 


J  jexp[-K(v,5',y)L]-exp 


-J  K(v,5(x),y(x))dx 


>dv  =  0 . 


(51) 


For  some  kinds  of  function  of  the  line  contour  instead  of  (51)  it  is  possible  to  require,  that  integrand 
expression  was  equalled  to  zero  at  any  wave  number,  namely 


exp 


[-k(v,S\y)l]: 


exp 


-J  K(v,5'(x),y(x))dx 


or 

k(v,5,  y )  =  j  J  k(v,SXx),  y(x))dx . 

Ll 

For  example,  for  the  Lorentzian  line  one  can  obtain  the  following  well  known  relation  (with  the  account  of 
Vj  □  max{y(x)}): 


S  y  =  —  J  S  (v)y  (x)dx .  (52) 

Ll 

The  situation  is  different  for  lines  with  exponential  contours.  It  is  not  here  possible  to  equate  to  zero 
integrand  function.  Therefore  for  averaging  half-width  it  is  necessary  to  use  an  additional  approximation. 


10-46 


RTO-EN-AVT-1 62 


NATO 

_ i _ 

OTAN 

Radiation  Modeling  in  Shock-Tubes  and  Entry  Flows 


The  integral  (51)  can  be  presented  in  the  following  form: 


2  _ 

J  \  exp  [-k(  v,  S ,  y )  Z-]  -  exp 


-  j  k(v,  S  (x),  y(x))dx 


>dv  =  0 , 


(53) 


where  v2  is  the  wave  number,  at  which  the  exponent  in  (53)  becomes  equal  to  zero  by  virtue  of 
exponential  reduction  of  absorption  coefficient  at  large  v.  Then,  started  from  the  theorem  about  average,  it 
is  possible  to  prove,  that  the  integrand  expression  will  be  equal  to  zero  though  in  one  point  v*  of  the 
spectral  range  [\’| ,  v2  ] ,  that  is 


k(v*,S,y)  =  jjK(v\S(x),Y(x))d* 


or 


5 

— exp 

2y 


f 


V 


* 

v 


if^W 

42y« 


r 

exp 

V 


*  'N 

V 


dx . 


Now,  one  can  replace  the  function  y(x)  in  the  exponent  by  its  average  value  y 


5.0  EXAMPLES  OF  THE  USE  OF  THE  RANDOM  MODELS  OF  MOLECULAR 
LINES  FOR  RADIATION  HEAT  TRANSFER  PROBLEMS  ARE 
PRESENTED  IN  SECTION  6.6.  5.0  THE  MONTE-CARLO  METHODS 


Monte-Carlo  simulation  methods  have  great  potential  for  creation  of  the  universal  computational  codes 
intended  for  prediction  of  spectral  emissivity  of  different  aerospace  objects.  First  of  all,  these  are:  light¬ 
scattering  plumes  of  rocket  motors  (signatures  of  rocket  plumes),  emissivities  of  entering  space  vehicles  in 
light-scattering  and  non-scattering  atmospheres,  meteors,  etc.  Performance  of  modern  computers  and 
parallel  computing  technologies  allow  solve  radically  the  main  problem  of  the  Monte-Carlo  imitating 
algorithms,  namely  reduce  dispersion  of  numerical  simulation  results. 


Basic  Monte-Carlo  imitative  algorithms  of  radiation  transfer  of  heat  radiation  in  arbitrary  inhomogeneous 
volumes  of  light-scattering  media  will  be  considered  in  sections  5. 1-5.8.  Next  a  peculiarity  of  the  Monte- 
Carlo  algorithms  to  solve  some  applied  problems  which  are  of  practical  interest  for  aerospace  applications 
will  be  considered  in  sections  5.9-5.10. 


Investigations  of  efficiency  of  different  simulation  algorithms  with  reference  to  problem  of  calculation  of 
spectral  signatures  were  performed  in  [25-28].  Analysis  of  five  simulation  algorithms  is  given  in  [25]: 

•  The  LBL-method  of  integration  of  radiation  heat  transfer  equation  on  spectrum  of  rotational  lines; 

•  Hybrid  statistical  model; 

•  Method  of  smoothed  coefficients; 
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•  The  two-group  method; 

•  The  LBL-method  with  small  number  of  trajectories  of  photons. 


Brief  description  of  the  algorithms  is  presented  below  (sections  5. 1-5.5).  Further  developments  of  the 
imitative  algorithms  are  presented  in  sections  5.6-5.10. 


5.1  Line-by-line  integration  of  radiation  heat  transfer  equation  on  spectrum 
of  rotational  lines 


The  given  spectral  range  is  divided  on  N  of  spectral  sub  regions  Acog.  In  limits  of  each  sub  region  the 
following  averaged  absorption  coefficient  is  introduced: 


-j-  j 

Aco„  r 


8  Aco  _  _ 


N, 


KP  (co)  +  KC  (co)  +  (co) 


dco, 


(1) 


where  k-  (co)  is  the  spectral  volumetric  absorption  coefficient  of  i- th  line;  N[  is  the  number  of  lines 
located  in  the  sub  region  Aco^  ;  kc  (co) ,  k_p  (co)  are  the  spectral  volumetric  absorption  coefficient  of  gas 
(in  continuum  spectra)  and  solid  (liquid)  particles. 


Let  us  assume,  that  all  lines  have  Lorentzian  contour,  therefore 


K i  (CO) 


L,i 


(co  — C00,;)  +YI,( 


(2) 


where  Shyu  are  the  intensity  and  half-width  of  i- th  Lorentzian  line;  co0?;  is  the  wavenumber  of  centre  of 
the  i- th  line. 

Imitative  Monte-Carlo  simulation  of  photon  trajectories  can  be  performed  for  each  spectral  group  Acog  as 
for  any  “grey”  medium  [28].  It  should  be  emphasize,  that  in  the  case  of  non-uniform  medium  coefficients 
Kg,  kp  ,  ,  k1  ,  Kj  ( co) ;  Si9  y  are  functions  of  coordinates. 


5.2  Hybrid  Statistical  Method  [26] 

The  idea  of  the  method  is  extremely  simple.  The  Monte-Carlo  imitative  algorithms  are  applied  not  for 
simulation  of  monochromatic  photon  groups,  but  for  simulation  of  the  propagation  of  photon  groups 
whose  energy  is  determined  by  averaging  over  spectral  band.  The  averaged  energy  of  photon  groups  for 
various  optical  paths  is  determined  with  the  use  of  random  models  of  real  linear  spectra.  The  Curtis  - 
Godson  method  is  used  to  take  into  account  inhomogeneity  of  the  optical  path.  Detailed  description  of  the 
random  models  is  presented  in  Chapter  3.  The  following  works  [17,19]  can  be  recommended  for  more 
deep  studying  of  the  models. 


This  approach  has  the  following  advantages: 

3)  The  integration  over  wave  number  in  this  approach  is  the  analytical  one,  and  it  can  be  performed  for  a 
separate  spectral  band  that  contains  from  one  to  hundreds  and  thousands  of  spectral  lines.  Another 
words,  it  means  that  the  spectral  problem  is  solved  just  once  for  the  entire  spectral  band. 

4)  To  implement  this  random  method  into  regular  Monte-Carlo  algorithm,  it  is  enough  to  modify  the 
subroutine  for  calculation  of  free  path  of  photon  groups  and  insignificantly  the  geometrical  module. 
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5)  The  type  of  the  spectral  lines  taking  into  account  (atomic  or  molecular),  their  location  in  spectral 
region  (regular  or  chaotic),  their  contours  (the  Lorenz,  the  Doppler,  the  Voigt)  are  of  no  significance. 
All  the  above-mentioned  peculiarities  are  taken  into  account  by  the  random  models  that  can  be  used 
without  any  restrictions. 

6)  If  it  is  necessary  to  take  into  account  strong  inhomogeneity  of  optical  paths,  the  Curtis  -  Godson 
approximation  can  be  applied  with  any  other  model  that  would  describe  the  environmental 
inhomogeneity  along  the  optical  path  more  adequately  [23,24]. 

7)  The  method  does  not  lose  any  advantages  as  compared  with  the  Monte-Carlo  simulation  method,  in 
particular,  in  respect  to  solution  of  problems  for  volumes  with  complex  geometries. 

The  drawback  of  the  hybrid  Monte-Carlo  simulation  algorithm  is  the  typical  for  the  classical  imitative 
Monte-Carlo  algorithms  and  random  models.  These  are: 

1)  High  computational  intensity  of  the  Monte-Carlo  algorithms; 

2)  The  intrinsic  error  of  random  band  models  which  does  not  allow  an  asymptotic  transition  to  precise 
line-by-line  integration. 

Let  us  now  consider  some  details  of  the  hybrid  method.  As  it  was  mentioned  above,  the  basis  of  the 
method  is  the  simplest  Monte-Carlo  simulation  method.  This  method  simulates  trajectories  of  photon 
groups  inside  each  homogeneous  element  of  volume,  which  is  introduced  by  finite  difference  grid 
overlying  the  whole  calculation  domain. 

We  consider  fixed  spectral  band  Acog  after  dividing  the  spectrum  under  consideration  into  such  bands. 
We  can  enter,  for  example,  ~  50  -f  250  cm'1  for  AcGg  in  the  infra-red  region  of  the  spectrum.  Spectral  lines 
of  molecules  and  continuous  (or  quasi-continuous)  absorption  spectra  may  exist  in  each  band. 

Calculation  by  the  hybrid  statistical  method  starts  from  the  formulation  of  the  position  vector  of  linear 
homogeneous  section  in  the  direction  of  motion  of  corresponding  photons  in  the  inhomogeneous  media 
under  consideration  (Z  -,  j  =  1,2, 3,..., N;  lx  =0;  lN  is  the  coordinate  of  the  point  of  intersection  between 

the  beam  and  the  boundary  surface).  The  nodes  of  the  calculation  grid  lj  correspond  to  the  points  of 

intersection  between  the  photon  direction  with  all  the  boundary  surfaces  of  homogeneous  volume 
elements. 

Each  homogeneous  volume  element  has  associated  coordinates  lj ,  vectors  of  gas  temperature  Tf , 
temperatures  of  particles  T? ,  molar  concentrations  of  optically  active  components  of  gas  mixture  . , 
concentrations  of  optically  active  liquid  and  solid  mixture  components  ■  (n%  j  is  the  number  density  of 

particles  of  given  dimensions,  if  the  mixture  is  dispersive,  then  ■  is  the  number  density  of  particles 
with  a  certain  average  dimension),  and  the  total  pressure  in  gas  mixture  /?  ■  .  The  dimensionality  of  vectors 
T? ,TP ,jc| ■  amounts  to  N- 1,  since  the  properties  inside  each  elementary  zone  are  assumed  to  be 
constant. 
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This  allows  calculations  of  the  following  optical  properties  that  determine  the  character  of  radiation 
transfer  along  the  beam  within  the  limits  of  each  section  of  vector  l- : 

3)  Coefficients  of  absorption  k v-  and  scattering  gJ  of  the  condensed  phase  (in  this  particular  case  the 

calculation  is  made  according  to  the  Mie  theory  [29]),  as  well  as  the  scattering  phase-  function  or 
average  scattering  cosine  ju;  that  characterizes  the  extent  of  asymmetry  of  the  phase  function.  The 

given  functions  are  considered  to  be  independent  of  wave  number  within  the  limits  of  spectral  band 
Aco g  ; 

4)  Absorption  coefficients  in  the  continuous  spectrum; 

5)  The  average  absorption  coefficient  in  spectral  lines  (S/d)kb  ; 

6)  Line  broadening  parameter  p^  • ; 

7)  Effective  pressure  (Pe  )k  h  .  . 

The  average  absorption  coefficient  in  the  spectral  lines,  the  spectral  line  broadening  parameter  and  the 
effective  pressure  are  calculated  for  each  molecular  band  b  of  each  component  k  that  contributes  to 
radiation  absorption  in  band  Acog  .  These  functions  were  calculated  using  the  wide  band  model  [22]. 

The  next  stage  of  the  algorithm  is  to  build  a  vector  of  effective  optical  thicknesses  for  each  band  of  the 
spectral  lines,  for  which  the  Curtis  -  Godson  method  is  used.  Let  us  assume  that  collision  broadening 
dominates  so  the  lines  have  a  Lorentzian  profile: 


where  p*.  is  the  density  of  the  k- th  component  of  the  mixture;  S ,  y  are  the  averaged  intensity  and  half¬ 
width  of  rotational  lines  in  the  b- th  lines  band  of  the  k- th  species. 


Assuming  that  the  principle  of  statistical  independence  of  transmission  in  spectral  line  bands  is  correct, 
and  accounting  for  continuum  absorption,  the  full  optical  thickness  is 


*5  = 


ZZ 

k  b 


Tk,bJ  + 


z 

k 


Tk,j  +, 


(6) 


10-50 


RTO-EN-AVT-1 62 


NATO 

i _ 

OTAN 

Radiation  Modeling  in  Shock-Tubes  and  Entry  Flows 


4,j  =  £<{L-L- 1 ) >  ^  =  £  Km  (Zm  - /«-! )  •  (7) 

m= 2  m=2 

If  the  Doppler  broadening  dominates,  the  calculation  of  Tk,b,j  is  changed  by  substituting  the  Doppler 
profile  for  Equation  (3)  (see,  for  example,  [19]). 

Simultaneously  with  calculation  of  the  optical  absorption  thickness,  calculation  of  optical  scattering 
thickness  is  performed  as  following 

T5  =  ZCTm(Z«-Zm- 1)>  (g) 

m-2 

which  enables  later  estimation  of  the  probability  of  absorption  (or  scattering)  of  photons  by  particles  of  a 
media. 

The  photons  mean  free  path  is  identified  as  usual.  However,  to  find  the  point  of  collision  between  the 
photons  and  particles  of  the  environment,  additional  calculations  have  to  be  made.  The  main  reason  for 
this  is  the  non-linear  relationship  for  transmission  in  spectral  groups  with  spectral  lines  (the  square -root 
low  for  the  equivalent  width  of  molecular  lines  [17]).  This  is  because  inequality  of  total  transmission  of 
the  entire  section  to  the  sum  of  individual  transmissions  of  its  components. 


Assume  that  after  generation  of  a  random  number  it  was  established  that  optical  thickness  xc  ,  that 
corresponds  to  the  mean  free  path,  satisfies  the  inequality: 


Ty  <TC  <Ty+1,  T*-  -Taj  +  T Sj. 


(9) 


It  is  then  necessary  to  find  coordinate  lc  corresponding  to  value  tc  .  In  a  general,  it  is  necessary  to 
organize  an  iterative  process  to  find  the  value  of  lc .  A  non-iterative  method  can  be  used  if  it  was  assumed 
that  Equation  (9)  corresponds  to  a  linear  relationship  between  the  optical  thickness  and  the  physical 
coordinate  which  is  accurate  for  low  optical  thicknesses  of  the  volume  elements.  Then: 


lc 


~lj  + 


(lj+i-lj)(xC-x'j) 


(10) 


Now  one  can  estimate  the  probability  of  absorption  (scattering)  in  the  point  where  the  photon  collides  with 
particles.  As  opposed  to  the  Monte-Carlo  calculation  of  monochromatic  radiation  transfer,  where  the 
probability  of  absorption  (scattering)  is  estimated  as  a  relation  of  the  absorption  (scattering)  coefficient  to 
the  full  attenuation  coefficient,  in  this  case  the  procedure  applies  to  corresponding  optical  thicknesses.  For 
instance,  the  probability  of  scattering  is  estimated  as  follows 


P.= 


lc-L 


rrS  I  n-a 

X  -f  T 

C  C 


Lj+ 1  J 


The  rest  of  the  Monte-Carlo  procedure  is  the  same  as  for  monochromatic  radiation. 


There  is  one  more  peculiarity  of  the  calculation  algorithm  with  molecular  lines.  It  is  connected  with  the 
problem  of  a  choice  of  energy  for  each  simulated  photon  group. 
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There  are  two  algorithms.  In  the  first  algorithm  the  energy  is  determined  as  uniform  for  all  photon  groups. 
This  algorithm  is  realized  in  the  following  steps: 

1.  The  integrated  emissivity  E  of  the  whole  researched  volume  V  and  separate  its  parts  En  (on  which 
the  volume  is  broken  by  the  finite  difference  grid)  are  calculated; 

2.  The  energy  of  each  photon  group  ef  is  determined  by  dividing  of  complete  energy  E  on  quantity 
of  simulated  groups  N  j- ; 

3.  The  relative  probability  pn  of  photon  emission  in  different  volume  elements  is  determined  by  the 
relation  of  the  emissivity  of  the  given  zone  to  integrated  emissivity:  pn=EnlE .  According  to 
this  relative  probability  number  of  photon  groups  emitted  by  the  given  volume  element  is 
determined  as  follows:  N f  n  =  pnN f  . 

So,  in  this  algorithm  all  photon  groups  have  identical  energy,  but  different  volume  zones  emit  various 
number  of  the  photon  groups. 

In  the  second  algorithm  each  photon  group  gains  different  energy,  determined  by  the  emissivity  of  given 
spatial  zone  (volume  element)  from  which  this  photon  group  is  emitted.  In  this  algorithm  it  is  expedient 
for  each  spatial  zone  to  determine  identical  number  of  emitted  photon  groups.  This  algorithm  is  realized  in 
the  following  steps: 

1 .  The  values  En  and  E  are  calculated; 

2.  The  number  of  photon  groups  emitted  by  each  spatial  zone  is  determined  under  the  formula 

Nf*=Nf/NV> 

where  Nv  is  the  number  of  spatial  zones. 

3.  The  energy  of  the  given  photon  group  is  calculated  under  the  formula 

ef  =  En/Nf,n  ■ 

Both  stated  algorithms  can  be  modified  by  weight  algorithm  of  modelling,  when  the  energy  of  groups 
photon  changes  during  collisions  [1]. 

The  mentioned  above  peculiarity  consists  in  a  way  of  estimation  of  energy  En  : 

E„=A-k\  j  k(co) J ^  gjdcodV7 ,  (11) 

Vn  Acog 

where  Vn  is  the  volume  of  n- th  spatial  zone;  Jb  is  the  spectral  intensity  of  the  Planck  radiation;  k(co)  is 
the  volumetric  spectral  absorption  coefficient.  As  properties  of  medium  are  averaged  in  limits  of  each 
zone,  and  the  spectral  ranges  Acog  are  usually  slender,  so  in  their  limits  with  good  accuracy  it  is  possible 
to  enter  average  spectral  Planck  intensity 


J b,n 


Aco 


J  J b,(n,n • 


8  Aco„ 


(12) 
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Then  the  expression  (1 1)  is  much  simplified 

En=^b,nVn  f  Kco ,«d®  •  (13) 

Aco? 

If  the  volumetric  spectral  absorption  coefficient  k„  does  not  depend  on  wave  number  in  limits  Acog 
(approximation  of  “grey”  medium),  then 


En  =  4nJb,n  ■  •  (14) 

When  spectral  region  AcOg  contains  TV/  spectral  lines,  then  integral  in  the  formula  (13)  may  be  calculated 
as  follows: 


Ni  Nt 

J  Kco,ndo)  =  K«Ac0g  +  KnAc0g  +  j  2X*  (“)  ~  K£A(0g  +  K«At0g  +2X«’  (15) 

Aco^  Aco^  i= 1  i= 1 

It  is  assumed  also,  that  the  half-width  of  lines  (with  spectral  volumetric  absorption  coefficient  k-  n  (co) ) 
satisfy  to  the  following  condition: 


Y;,JP  Acog, 


then 


Si,n  =fKL(®)dc0~  j  Kln(®)d®- 


(16) 


0  Acog 

Thus,  the  integrated  emissivity  is  calculated  as  follows: 


En  =^Jb,nVn 


N,  \ 

i=l  J 


(17) 


Formulas  (13)  and  (17)  are  easily  generalized  on  the  case,  when  the  temperatures  of  gas  and  solid  (liquid) 
phases  differ: 


En=4KVn(jbpnKZ  +  jenKs)Ao>g, 


or 


E„  =  4nV„ 


JbPn<A(°g 


+  J8 

±Jb,n 


N, 


<A®g  +  2X« 


(18) 


where  j£n ,  Jjjn  are  the  Planck  intensities  at  temperatures  of  particles  Tp  and  gas  T8  ;  kp  ,  k|  are  the 
overage  volume  absorption  coefficients  of  particles  and  gas. 


It  is  necessary  to  pay  attention  on  the  fact,  that  the  stated  way  of  estimation  of  emissivity  is  well 
reasonable  in  cases,  when  the  volume  Vn  is  optically  thin.  If  the  optical  thickness  xn  of  the  volume  Vn 
does  not  satisfy  to  this  condition,  i.e.  xn  >  1 ,  the  specified  way  of  calculation  of  emissivity  can  result  in 
significant  errors  (in  particular  by  use  of  a  hybrid  Monte-Carlo  method). 
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5.3  Method  of  the  smoothed  coefficients 

This  method  is  actually  equivalent  to  the  hybrid  Monte-Carlo  method,  but  here  the  following 
simplification  is  used 


(19) 


^ k,b,j  %k,b,j  * 

This  relation  is  correct  under  the  following  condition: 


It  is  obvious,  that  with  increase  of  optical  thickness  it  is  necessary  to  expect  increase  of  the  error  of  the 
method.  However,  as  it  was  shown  in  [25],  at  very  large  optical  thickness  this  method  gives  more 
physically  reasonable  result,  than  hybrid  Monte-Carlo  method. 

5.4  The  two-group  method 

Let  spectral  range  AcOg  contains  TV/  molecular  lines.  Let  us  assume,  that  integrated  emission  in  the  given 
spectral  range  does  not  depend  on  location  of  the  molecular  lines. 

Average  intensity  and  half-width  for  all  lines  from  the  range  Acog  in  each  spatial  zone  is  calculated  as 


follows: 


(20) 


Then  the  average  half- width  on  all  spatial  zones  is  determined  as  follows: 


(21) 


Let  us  enter  two  spectral  ranges  inside  the  spectral  range  Acog  such  that: 


Acc*!  =  min  j  2 Nl  (y ) ,  Acog  J , 

Aco2  =  Aco^  -  Aooj. 


(22) 


The  average  absorption  coefficient  in  the  spectral  ranges  Acc^  and  Aco2  are  determined  under  the  formulas 


(*n) 


^2 ,n  ^ n 


(23) 


7i(y)’ 


Now  everything  is  ready  to  imitative  simulation  of  the  two  group  photons  with  the  averaged  absorption 
coefficients  (23).  Emissivity  of  each  volume  elements  is  determined  in  this  case  as  follows: 

•  En=  4nVnJb ^k^Aooj  -  for  the  first  spectral  region; 

•  En=  4nVnJh ?nk2  n Acd2  -  for  the  second  spectral  region. 
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5.5  Line-by-line  integration  with  little  number  of  trajectories 

Estimations  of  dispersion  of  direct  statistical  simulation  results,  and  also  experience,  accumulated  at  the 
solution  of  similar  problems  in  “grey”  statement  show,  that  for  obtaining  of  satisfactory  accuracy  of  the 

simulation  it  is  enough,  as  a  rule,  to  simulate  about  Nf  =104  -^5xl04  trajectories.  For  the  case  of  relative 
low  albedo  ( co  =  — - —  <  0.9 )  with  moderate  optical  thickness  ( x  <  1 )  this  number  can  be  reduced  yet  by 

CT  +  K 

the  order.  And  on  the  contrary,  for  multi-scattering  medium  of  large  optical  thickness  this  number  is 
necessary  yet  to  increase. 

When  there  is  a  need  to  find  averaged  over  the  spectral  region  AcOg  radiating  characteristics,  it  is  possible 
to  use  the  following  algorithm.  Let  us  try  to  find  the  solution  of  the  problem,  simulating  trajectories  only 
of  Nf  ~  104  photon  groups.  The  energy  of  each  new  photon  group  will  be  estimated  statistically.  For  this 
purpose  spectral  absorption  coefficient  is  calculated  in  any  spectral  point  inside  the  region  AcOg  . 

The  stated  algorithm  can  be  modified  as  follows.  We  shall  break  the  range  AcOg  on  elementary  spectral 
sub  regions.  The  number  of  these  sub  regions  must  be  sufficient  for  the  detailed  description  of  line 
structure  of  the  spectrum.  For  example,  at  Aa>g  =20  cm1  and  average  size  of  half-width  y  =  0.1  cm1  it  is 
enough  to  enter  ~  1000  spectral  sub  regions  AcOg  . 

Further  to  calculate  averaged  radiating  characteristics  one  can  apply  usual  procedures  of  line-by-line 
integration,  but  within  the  bounds  of  the  each  spectral  sub  region  AcOg  to  simulate  trajectories  not  of  ~  104 
photon  groups  (as  in  the  regular  line-by-line  method),  but  only  -10-15  photon  groups.  It  should  be 
stressed,  that  it  is  impossible  to  determine  the  spectral  characteristics  inside  AcOg  by  this  way,  but  the 
averaged  characteristics  in  the  region  AcOg  can  be  obtained  with  rather  good  accuracy. 

5.6  The  Monte-Carlo  imitative  method  based  on  the  Maximum  Cross  Section 
(MCS)  method 

Regular  imitative  Monte-Carlo  algorithms  considered  above  were  developed  for  calculation  of  averaged 
radiation  transfer  over  rotational  lines  in  spectral  range  AcOg  on  the  orthogonal  calculation  grids  in  2D 
cylindrical  geometry.  These  imitative  algorithms  were  used  for  test  calculations  of  signatures  of  model 
rocket  motor  plumes  [28]  and  for  comparison  with  computational  data  [30] .  These  algorithms  were  used 
also  at  study  of  different  optical  models  for  prediction  of  rocket  motor  plume  signatures  [27].  It  was 
shown,  that  considered  imitative  Monte-Carlo  algorithms  can  be  used  quite  really  even  with  “line-by-line” 
(LBL)  models  of  rotational  line  structure.  Development  of  these  imitative  Monte-Carlo  algorithms  with 
reference  to  non-orthogonal  grids,  which  are  usually  used  at  the  solution  of  gasdynamic  problems,  is 
presented  in  the  Ref.  [31].  Three  imitative  algorithms  on  non-orthogonal  grids  were  considered  there: 

8)  The  LBL-method  based  on  the  imitative  Monte-Carlo  algorithm  and  the  MCS  method; 

9)  The  LBL-method  based  on  the  imitative  Monte-Carlo  algorithm  and  the  method  of  quasi-random 
sampling  for  formation  of  straight-line  trajectory  segments  (further  we  will  use  also  the  term  “direct 
Monte-Carlo  simulation”); 

10)  Hybrid  method  based  on  the  method  of  quasi -random  sampling. 
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The  Maximum  Cross-Section  method  is  based  on  the  obvious  statement  that  adding  any  quantity  of  8- 
scatterer  (8-diffuser  of  light)  into  investigated  volume  of  the  light-scattering  media  does  not  change  the 
process  of  photons  diffusion  (by  the  definition,  the  S-scatterer  is  the  scatterer  without  changing  direction). 

Let  us  introduce  the  maximum  total  extinction  coefficient  of  the  multi -component  multiphase  polydisperse 
medium,  using  the  following  definition: 

Pr,max  —  iflUK  jp^  1 1 ,  P t,l  ~  ^ g,l  ^p,l  ’  (24) 

where  l  is  the  number  of  the  volume  elements  defined  by  four  nodes  indexes  of  introduced  calculation 
grid;  are  the  gas  and  particles  absorption  coefficients;  g/  is  the  averaged  scattering  coefficient 

of  polydisperse  solid  phase;  $tj  is  the  extinction  coefficient. 

For  each  /- th  elementary  volume  of  the  computational  grid  there  may  be  introduced  also  the  8-scattering 
coefficient: 


P5 ,/  -Pr,max  “P t,l  •  (25) 

The  use  of  introduced  functions  ( Pr?max  >  P t,i  >  and  P§,/ )  allows  simulate  the  photons  diffusion  in  a 
heterogeneous  medium  as  well  as  in  a  homogeneous  one.  In  so  doing  it  is  agreed  that  at  each  photon’s 
collision  with  particles  of  the  medium  the  8-scattering  occurs  with  the  probability  of  P§  /  /P^max  >  and  the 
real  interaction  occurs  with  the  probability  of  1-P§,//P?,max  •  Another  words,  in  the  maximum  cross- 
section  method  the  delta-scattering  probability  is  determined  as  =P5/Pr,max  >  and  the  real  interaction 
probability  is  determined  as  pS  K  =  1  -  p§ .  Thus  for  simulation  of  collision  type  one  should  check  the  fact 
of  satisfaction  of  the  following  condition  pd>  y  ( y  is  the  quasi-random  number  from  the  equally 
probable  distribution  in  segment  [0,1]).  In  this  case  the  delta-scattering  has  been  occurred.  Otherwise  the 
true  collision  has  been  occurred. 


The  true  collision  may  be  both  absorption  and  scatter.  Generating  new  random  number  y  and  comparing 
it  with  the  scattering  probability  ps  =  a/pf  one  can  make  a  conclusion  concerning  the  type  of  the 
collision.  If  y  <  ps  there  has  been  occurred  scattering;  and  if  y  >  ps  or,  what  is  the  same,  y  <  pK  =  1  -  ps 
it  is  considered  that  there  has  been  occurred  absorption. 


The  use  of  the  maximum  cross-section  method  with  the  transport  approximation  of  scattering  processes 
simplifies  simulation  algorithm  to  a  greater  degree.  The  total  “delta-scattering”  coefficient  can  be 
presented  in  these  cases  as  follows: 

p6=P/,max-[K  +  (1-pH  Or  PS  =  Pt,max  “  P?  +  • 


Another  word,  in  the  transport  approximation  anisotropy  of  any  scattering  may  be  taken  into  account  at 
the  level  of  collision  type  modeling  (true  or  not  the  given  collision).  If  the  following  inequality 


P8,s 


g-g(I-li)  _ 
- - — —  =  P>Y, 

G 


is  true,  then  this  is  the  ^scattering;  and  if  p6  s  <  y  is  true,  then  this  is  the  real  scattering. 


5.7  Imitative  Monte-Carlo  algorithm  based  on  the  quasi-random  sampling  of  photon 
trajectory  parameters 

To  calculate  radiation  transfer  through  any  optical  path  with  directional  vector  Cl  there  is  necessity  to 
know  optical  properties  along  this  path.  For  this  purposes  the  coordinates  of  the  ray  Cl  intersections  with 
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all  surfaces  of  the  spatial  finite-difference  mesh  should  be  calculated.  This  is  illustrated  in  Figure  3.9. 
Relations  of  the  analytical  geometry  should  be  used  for  these  purposes.  As  it  was  shown  above,  at 
consideration  the  ray-tracing  method  (see  section  3.3),  this  algorithm  is  not  effective  due  to  very  high 
laboriousness,  especially  for  non-structured  meshes.  Algorithm  of  the  quasi-random  sampling  can  be  used 
with  the  Monte-Carlo  imitative  algorithms,  as  it  was  used  with  ray -tracing  method. 

5.8  The  hybrid  method  based  on  the  quasi-random  sampling  method 

This  is  the  hybrid  simulation  method  (see  section  5.2  and  [25,26]),  which  is  united  in  the  given  case 
together  with  the  method  of  quasi-random  sampling. 

5.9  Three-dimensional  simulation  algorithms 

Imitative  algorithms  considered  above  can  be  applied  for  prediction  of  spectral  emissivity  of  axially- 
symmetric  volumes,  for  examples,  rocket  plumes.  Calculation  of  the  directed  radiating  ability  of  a  single 
plume  is  of  significant  practical  interest  for  comparison  of  the  results  with  experimental  and  calculation 
data  of  other  authors.  The  overwhelming  majority  of  all  calculations  were  based  on  the  model  of  axially 
symmetric  plume  [26-28,31].  But  very  often  radiating  aerospace  objects  must  be  considered  as  the  three- 
dimensional  objects.  For  example,  heat  radiation  of  single  and  multi -block  plumes  or  light  scattering 
volumes  illuminated  by  sources  of  external  radiation,  for  example  sunny  [32].  It  is  presumed  in  this  case, 
that  the  plane-parallel  flux  of  radiation  illuminates  one  of  the  planes  of  calculation  area  under  any  angle  of 
attack.  Heat  radiation  of  Sun  can  be  compared  to  radiation  of  absolutely  black  body  at  temperature 
~  6000  K.  Spectral  dependence  of  scattered  sunlight  is  formed  under  action  of  absorption  of  gas  and  solid 
phases  of  the  plume. 

Spectral  radiating  ability  of  multi -block  jets  even  in  relation  to  self -radiation  is  the  three-dimensional 
problem.  Three-dimensional  imitative  Monte-Carlo  codes  were  developed  and  tested  in  [33].  Figures 
6.40-6.44  of  section  6.7  show  example  of  distribution  of  gas  temperatures  and  concentration  of  condensed 
phase  in  three-block  plumes.  Examples  of  numerical  simulation  results  for  such  a  configuration  of  multi - 
block  radiating  plume  are  also  presented  in  the  section  6.7. 

5.10  Monte-Carlo  Local  Estimation  of  Directional  Emissivity  (MCLEDE) 

The  Monte-Carlo  method  considered  above  is  useful  when  it  is  necessary  to  find  the  angular  distribution 
of  radiation  intensity  over  whole  solid  angle.  In  this  case  each  photon  leaving  the  radiating  volume  adds  its 
contribution  into  the  number  of  the  photons  moving  in  a  concrete  zone  of  calculation  grid  over  angular 
variables,  and  at  the  contraction  of  this  angular  zone  it  is  necessary  to  increase  the  number  of  modelled 
photons  so  that  the  statistical  error  of  the  calculation  won't  increase.  However,  in  practice  it  is  often 
necessary  to  know  the  radiation  intensity  in  the  strictly  determined  direction.  In  that  case  it  is  possible  to 
speed  up  the  performance  of  the  calculation  code  considerably,  at  the  same  time  having  made  the  angular 
range  infinitesimal.  Numerical  simulation  algorithm  for  such  a  case  was  presented  in  [34].  Here  we  will 
consider  a  substantiation  of  the  algorithm. 

Suppose  we  need  to  find  radiation  intensity  in  any  strictly  determined  direction  \ix .  We  shall  establish 
counters  not  for  the  number  N  of  the  photons  moving  in  the  range  dp ,  but  for  their  density  A/dp .  For 
each  photon  we  shall  set  initial  energy  E0.  After  each  collision  the  energy  of  a  photon  will  decrease: 
E  =  coEq.  That  corresponds  to  the  estimation  of  probability  for  the  given  photon  to  continue  the  movement 
in  any  direction  (for  the  case  of  isotropic  scattering). 
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If  we  count  the  photons  moving  in  the  range  [pj  -  dp;  \ii  +  dp] ,  the  photon  with  energy  £dp  will  move 
in  the  direction  p! ,  and  the  energy  of  the  photons  moving  in  all  directions  p2  ^  Pi  will  be  ^-(1-dp) . 
The  direction  p2  is  determined  by  any  algorithm  of  scattering  simulation.  After  that  it  is  necessary  to  trace 
both  of  these  photons. 

The  photon  with  the  energy  £dp  =  mE,0dp  moves  in  the  direction  p}  before  the  following  collision 
occurs.  After  this  collision  the  energy  of  the  photon,  moving  along  \ix ,  becomes  equal  to  oo2^  (dp)2  .  All 
the  counters  in  computational  cells  between  these  collisions  increase  by  coE0 .  After  the  next  collision  the 
counters  will  increase  by  co^odp ,  and  at  dp  — »  0  we  can  neglect  this  increase.  Due  to  this  fact,  we  shall 
trace  the  photon  with  energy  Ed\i  only  until  the  first  collision  occurs.  At  the  same  time  we  shall  trace  the 
photon  with  the  energy  £(l-dp)  (moving  in  the  direction  p2)  by  a  common  algorithm.  This  photon 
won't  affect  the  counter  directly,  but  will  create  new  photons  with  energy  £(l-dp)dp  «  £dp  after  each 
collision.  The  radiation  intensity  corresponding  to  the  value  of  the  counter  E%  in  the  given  cell  is 
/  =  EZ/ / 4ny.i  . 

It  should  be  stressed  that  if  we  want  to  find  radiation  intensity  only  on  the  layer  border,  we  don't  need  to 
determine  the  free  path  of  the  photon  with  the  energy  £dp  .  Instead,  we  can  just  multiply  this  energy  by 

Let's  estimate  the  ratio  of  the  performance  of  the  probabilistic  algorithm  t\  to  the  performance  of  the 
regular  Monte  Carlo  algorithm  (at  the  statistical  error  8  fixed).  It  is  obvious  that 


hlh  ~  N2/Nl ,  (21) 

where  N\  and  N2  are  the  numbers  of  the  photons  required  by  the  probabilistic  and  the  common  methods 
correspondingly  to  reach  the  error  e  .  This  error  can  be  given  by 

z  =  c/4n  ,  (22) 

where  n  is  the  number  of  the  photons  which  have  made  the  contribution  to  our  statistics,  and  c  is  a 
constant,  almost  identical  for  both  algorithms  at  small  dp  .  For  the  probabilistic  method  rt\  =  N\ ,  and  for 
the  common  algorithm  n2  =Af2dp-/(p),  where  /(p)  is  a  factor  with  the  meaning  of  about  1.  Having 
assumed  /(p) « 1 ,  from  (21)  we  shall  have 


hlh  ~  1/dp-  (23) 

In  tree-dimensional  case  we  shall  have  (24)  instead  of  (23): 

t2/t1  ~  1/  dpdcp .  (24) 

To  compare  regular  algorithm  with  the  local  estimation  algorithm  a  special  test  case  was  considered  [34]. 
The  specified  test  case  has  been  offered  in  work  [30].  The  radiating  cylinder  is  filled  with  a  homogeneous 
high-temperature  mixture  of  gases  H20,  C02,  CO  at  temperature  T  =  2000  K  and  with  particles  A1203  with 
an  average  radius  of  2  microns.  The  radius  of  the  cylinder  is  taken  equal  to  10  cm,  and  its  height  to  600 
cm.  For  calculations  of  spectral  optical  properties  the  optical  model  [35-39]  of  the  molecular  lines 
averaged  on  rotary  structure  was  used. 

Figure  5.1  shows  the  spectral  radiation  signature  of  the  cylinder  without  scattering  at  the  observation  angle 
of  90  degrees.  Calculations  were  performed  using  two  imitating  algorithms. 
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Figure  5.1 :  Spectral  signatures  predicted  by  two  different  Monte-Carlo  algorithms.  No  scattering.  Results 
obtained  by  the  regular  Monte-Carlo  algorithm  are  presented  for  the  two  angles  of  observation 

One  of  them  the  Monte-Carlo  local  estimation  of  directional  emissivity  algorithm  and  another  one  is  the 
regular  Monte-Carlo  imitative  algorithm  [26,28].  The  method  of  the  local  estimation  of  emissivity  has 
shown  high  efficiency.  This  method  reaches  5%  accuracy  of  statistical  modeling  approximately  100  times 
faster. 

Figures  5.2  show  spectral  directional  emissivity  at  solid  particle  concentration  107cm-3.  In  this  case, 
histories  of  104  photons  were  simulated.  One  can  see  advantage  of  the  local  estimation  method  for  the 
problem  under  consideration.  Details  of  these  calculations  are  presented  in  [34]. 


Spectral  directional  emissivity,  W/(sr*mcm)  I  Plane  X~z] 


Figure  5.2:  Spectral  signatures  of  light-scattering  volume  predicted  by  two  different  Monte-Carlo  algorithms 
at  number  of  photon  histories  104 .  Results  obtained  by  the  regular  Monte-Carlo  algorithm  are  presented  for 

the  two  angles  of  observation 
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6.0  EXAMPLES  OF  APPLICATION  OF  METHODS  FOR  SOLVING 
RADIATION  HEAT  TRANSFER  EQUATIONS 

6.1  Pi-approximation  [4] 

The  radiation  heat  transfer  problem  was  considered  for  two-dimensional  axially-symmetric  geometry  to 
analyze  typical  solutions  of  the  Prequations.  Temperature  distribution  in  such  a  volume  is  shown  in 
Figure  6.1.  Note  that  such  a  spherically  symmetric  radiating  region,  localized  in  space,  allows  test  wide 
class  of  methods  developed  in  the  radiation  heat  transfer  theory. 
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Figure  6.1 :  Temperature  inside  two-dimensional  cylindrical  volume 


Temperature  field  shown  in  Figure  6.1  was  calculated  by 


T{X’r)  =  Tmin+{Trmx 


f 

Tmm  )exP 

V 


Rm] 


where  7^  =  300  K,  I ^  =  1 8000 K,  m-  4, 

7?  =  ^/(x-^)2+(r-r0)2  , 
v()  =5cm,  r0  =  0,  R()  =  1  cm  . 

Functions  of  radiation  field  corresponded  to  this  temperature  distribution  and  k  =  0. 1  cm'1  obtained  by  the 
Pi-approximation  (see  Equation  (3.40))  are  shown  in  Figures  6. 2-6.5. 
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r,  cm 


Figure  6.2:  Volume  density  of  radiation,  J/cm3;  K  =  0.1cm  1 
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Figure  6.3:  Radiation  flux  Wx,W/cm2;  K  =  0.1cm  1 
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Figure  6.5:  Divergency  of  radiation  flux  Qr  =div  W ,  W/cm2;  K  =  0.1cm  1 


Finite-difference  grid  used  at  these  calculations  is  shown  in  Figure  6.6.  Calculated  data  presented  here 
may  be  used  at  verification  of  other  methods  and  codes  of  radiation  heat  transfer  theory. 
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Figure  6.6:  Finite  different  grid 


6.2  The  quadro-moment  method  [40] 

The  quadro-moment  method  was  applied  to  integration  of  two-dimensional  radiation  heat  transfer  problem 
similar  to  those  considered  in  the  previous  section.  Temperature  distribution  with  the  centre  located  at  the 
symmetry  axis  of  cylindrical  volume  was  calculated  as  follows  (Figure  6.7): 


T  (  x,  r)  =  T0  +  Tx  •  exp 


r  \2 


\hrj 


+ 


x-x,. 


\2 


bx  J 


where  T0  =300K,  7]  =10000K ;  xc  is  the  axial  coordinate  of  the  high  temperature  area  centre,  br,bx  are 
the  shape  form  parameters  of  the  temperature  distribution  along  axes  r  and  x. 
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Figure  6.7:  Temperature  distribution 
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The  radius  and  height  of  the  cylinder  were  entered  as  5  cm.  These  calculations  were  performed  with  the  Pi 
approximation  of  the  spherical  harmonics  method  and  with  the  zeroth-order  approximation  of  the  quadro- 
moment  method.  Different  calculation  grids,  absorption  coefficients  k  and  coefficients  defining  the  value 
of  artificial  calculation  diffusion  s  were  used.  At  all  figures  representing  the  isometric  lines  of  divergence 
of  the  radiant  flux  integral  at  0.5  -fO.6  pm,  the  figures  at  the  curves  denote  the  values  of  Q  in  W/cm3. 
Symbol  (+)  here  means  emission  of  radiation,  while  the  symbol  (-)  means  absorption. 


Comparison  of  calculated  data  presented  in  Figure  6.8.  It  is  clear  from  this  figure  and  it  was  shown  [40] 
that  the  zero  approximation  of  the  quadromoment  method  allows  obtain  results  close  to  those  obtained 
with  the  Pi -approximation  of  the  spherical  harmonics.  But  it  should  be  taken  into  account  that  the  quadro¬ 
moment  method  has  some  advantages  as  compared  to  Prapproximation  at  the  use  of  the  methods  with 
multi-group  spectral  method. 
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Figure  6.8:  Radiation  flux  divergence:  Pi-method  (left),  the  quadro-moment  method  (right);  K  =  0.1cm_1 
6.3  The  ray-tracing  method 

6.3.1  Radiative  heating  of  internal  surface  of  the  laser  plasma  generator  [41,42] 

Radiation  heat  transfer  was  considered  in  working  zone  of  laser  plasma  generator  (LSPG).  Calculations 
were  performed  for  the  following  input  data:  the  capacity  of  the  CW  C02-laser  is  PL  =  200  kW,  the  focal 
length  of  the  focusing  lens  is  /  =  3  cm,  the  divergence  of  the  laser  radiation  is  0  =  0. 1  mradian,  the  initial 
radius  of  the  laser  beam  is  Rb  =  1.0  cm,  the  pressure  in  the  LSPG  channel  is  p0  =  1  atm,  the  gas  velocities 
in  the  input  section  is  u0=  30  m/s,  the  length  and  radius  of  the  LSPG  channel  are  Lc  =  11cm  and 
Rc  =1.3  cm.  Calculation  domain  was  covered  by  inhomogeneous  calculation  grid  NJ  =120,  NI  =  40  (NJ 
is  the  number  of  calculation  grid  nodes  along  %-axis;  NI  is  the  number  of  nodes  along  r-axis). 

Temperature  profile  is  shown  in  Figure  6.9  for  laminar  gas  flows  in  LSPG  [42]. 
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Figure  6.9:  Temperature  distribution  (in  K)  at  entrance  velocity  u0  =30  m/s  and  PL  =  200  kW  in  air  LSPG; 

the  laminar  gas  flow 


Integral  radiation  flux  presented  in  Figure  6.10  calculated  for  148-group  spectral  model.  The  ray-tracing 
method  was  used  for  prediction  of  the  group  spectral  fluxes  on  the  internal  surface.  Example  of  such 
calculations  is  presented  in  Figure  6.11,  where  the  148-group  radiation  fluxes  are  presented  for  the  five 
points  on  the  LSPG  internal  surface.  Numbers  of  presented  points  correspond  to  the  following  coordinates 
along  internal  surface:  1  -  x  =  0.46  cm,  2  -  x  =  0.93  cm,  3  -  x  =  1 .79  cm,  6  -  x  =  4. 1  cm,  11  - 
x  =  10.8  cm. 


Integral  radiative  flux,  W/cm**2 


X,  cm 

Figure  6.10:  Integral  radiation  heat  flux  along  internal  cylindrical  surface  of  hydrogen  LSPG,  W/cm2. 
The  ray-tracing  method.  Number  of  angular  points  41.  Number  of  spectral  groups  148 
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Figure  6.11 :  Group  radiation  heat  fluxes  at  different  points  of  internal  surface  of  air  LSPG,  W/cm2. 

The  Ray-tracing  method.  Number  of  angular  points  41 .  Number  of  spectral  groups  148 

6.3.2  Radiative  heating  of  space  vehicle  surface  [43,44] 

Computational  code  NERAT-2D  was  used  for  prediction  of  radiative  heating  of  whole  surface  and 
spectral  signature  of  space  vehicles,  entering  into  planet  atmospheres.  The  models  which  are  realized  in 
the  computational  code  NERAT-2D  include  description  of  general  physical-chemical  processes  which  can 
be  significant  for  solving  problem  under  consideration,  namely:  the  complicated  gasdynamic  structure 
appearing  at  streamline  of  space  vehicles  of  complex  geometry  is  described  by  the  Navier-  Stokes 
equations  for  viscous  heat  conducting  gas;  the  non-equilibrium  physical  and  chemical  kinetics  of  gas 
behind  shock  wave  and  in  wake  are  described  by  the  equations  of  energy  conservations  for  different 
internal  degrees  of  freedom  and  by  the  species  mass  conservation  equations;  the  spectral  radiation  transfer 
in  whole  disturbed  area  is  described  by  the  Prapproximation  and  by  the  ray-tracing  method. 

Prediction  of  the  radiation  heating  of  MSRO  surface  was  performed  in  [43,44]  for  four  trajectory  points, 
and  for  two  assumptions  concerning  the  space  vehicle  wall  catalicity  effects.  In  the  first  case  the  surface 
was  presumed  as  non-catalytic,  and  the  pseudo-catalytic  surface  with  97%  C02  and  3%  N2  at  the  wall  was 
presumed  for  the  second  case. 

Examples  of  numerical  simulation  results  are  presented  in  Figures  6.12-6.14.  These  calculation  data 
correspond  to  the  following  trajectory  point  in  Martian  atmosphere  (see  Table  6.1  and  [43,44]): 
Poo  =3.07  xlO_7g/cm3,  Poo  =82.3  erg/cm3,  Vo0  =3998  m/s.  Total  radiation  flux  on  MSRO  surface  is 
shown  in  Figure  6.13.  Figure  6.13  shows  also  the  spectral  radiation  fluxes  at  five  points  on  the  surface. 
Location  of  marked  points  one  can  determine  by  Figure  6.14.  Also  radial  and  axial  coordinates  of  marked 
points  (in  cm)  are  shown  on  the  field  of  each  figure.  Legends  for  the  figures  showing  total  and  spectral 
radiation  fluxes  contain  also  the  following  data:  Ng  is  the  number  of  spectral  groups;  Nm  is  the  number 
of  points  along  each  ray  emitting  from  MSRO  point;  Nq  is  the  minimal  number  of  rays  emitting  in  the  0- 
angle  direction  (0  is  the  angle  of  latitude  in  local  coordinate  system  connected  with  an  elemental  area  on 
MSRO  surface);  N 9  is  the  minimal  number  of  rays  emitting  in  the  cp-angle  direction  (cp  is  the  azimuth 
angle  in  local  coordinate  system).  It  should  be  noted  that  a  special  numerical  study  was  performed  for 
substantiation  of  the  chosen  numbers  of  angle  directions. 
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Figure  6.12:  Temperature  distribution  (left)  and  mass  fraction  of  CO2  (right).  Trajectory  point  No.3. 

Pseudo-catalytic  surface 
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Figure  6.13:  Total  radiation  flux  on  the  MSRO  surface  (left)  and  the  spectral  radiation  fluxes  at  some  points  on 
the  MSRO  surface.  Trajectory  point  No.3.  Pseudo-catalytic  surface.  Ng=9\  A/m  =  100,  A/e=11,  A/(p  =  11 


Figure  6.14:  Location  of  marked  points  on  the  MSRO  surface,  and  correspondence  of  point  numbers  to  axial 

coordinates  on  the  MSRO  surface 
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6.3.3  Spectral  signatures  of  space  vehicle  MSRO  [45] 

Spectral  signatures  of  heated  gas  behind  shock  wave  and  from  wake  at  MSRO  entering  into  Martian 
atmosphere  were  calculated  for  two  angles  of  observation  (see  schematic  of  the  problem  in  Fig.  6.15), 
0  =  30°  and  90°  .  Calculations  were  performed  for  four  trajectory  points.  The  ray-tracing  method  was  used 
in  the  case. 


Table  6.1 :  MSRO  trajectory  points 


No. 

of  trajectory 
point 

Time,  s 

Poo ,  g/cm3 

Poo  9 

erg/cm3 

Voo,  m/s 

Tm,  K 

i 

70 

3.14xl0-8 

8.4 

5687 

140 

2 

115 

2.93  xlO"7 

78.7 

5223 

140 

3 

175 

3.07  x  10-7 

82.3 

3998 

140 

4 

270 

2.82  xlO-8 

7.6 

3536 

140 

Presented  data  allow  confirm  well  known  fact  that  that  most  intensive  radiation  in  spectral  signature  is 
located  in  the  UV  spectral  region  for  high-velocity  trajectory  points  (Figure  6.16,  a)  and  in  the  infrared 
spectral  region  for  relatively  low-velocity  trajectory  points  (Figure  6.16,  d).  Infrared  radiation  generally 
connected,  as  before,  with  vibrational  bands  of  C02  and  CO. 

It  should  be  noted  that  ablation  of  heat  shields  of  lander  may  make  a  significant  contribution  to  the 
radiation  of  the  shock  layer  and  especially  of  a  wake.  There  are  two  consequences  from  this  fact.  First  of 
all,  the  model  must  include  detailed  description  of  the  thermo  chemical  processes  on  a  streamlined  surface 
and  predict  chemical  composition  and  mass  of  products  of  heat  protection  material  thermo  chemical 
destruction.  Secondly,  this  fact  can  be  used  for  increase  of  the  spectral  signature  at  any  given  spectral 
region  by  including  special  luminous  elements  as  admixture  of  heat  protection  material. 


Figure  6.15:  Temperature  field  for  trajectory  point  No.1 
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Figure  6.16:  Wavenumber  dependence  of  total-in-groups  signature  predicted  for  0  =  30°  . 

Trajectory  points  No.1  (a),2  (b),3  (c),4  (d).  The  pseudo-catalytic  surface. 

6.4  Discrete  ordinates  method  [16] 

Discreet  ordinates  method,  described  in  Section  3.4,  was  used  for  after  body  radiative  heating  prediction 
for  MSRO  space  vehicle.  It  was  assumed  that  Martian  atmosphere  contents  97%  C02  and  3%  N2.  The 
optical  region  was  chosen  as  1970-^4000  cm-1.  The  temperature  distribution  (Figure  6.17)  and  mass 
fraction  behind  space  vehicle  were  taken  from  [43,  46].  Surface  of  the  MSRO  was  assumed  absolutely 
black  at  the  temperature  500  K.  To  define  group  absorption  coefficients  the  ASTEROID  code  was  used  [1, 
47].  For  verification  of  the  radiation  heat  flux  prediction  the  NASA  standard  infrared  radiation  model  [35] 
was  also  used.  Coordinates  of  the  cylindrical  cell  are  defined  by  two  spatial  coordinates  r,  h  and  one 
angular  coordinate  0  .  Spatial  mesh  for  0  =  0  is  shown  in  Figure  6.18.  The  tetrahedral  grid  was  generated 
by  subdividing  cylindrical  cells  into  tetrahedrons.  Integral  radiation  heat  fluxes  along  afterbody  surface  of 
the  MSRO  depending  on  order  of  used  quadrature  (SN)  and  on  number  of  spectral  groups  (GR)  are  shown 
in  Figures  6.19  and  6.20.  Analysis  of  these  and  others  numerical  simulation  results  shows  that  developed 
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in  [16]  version  of  the  discreet  ordinates  method  allows  predict  radiative  heating  of  space  vehicle  surface 
with  good  accuracy. 
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Figure  6.17:  Temperature  distribution  behind  the  MSRO 


Figure  6.18:  Computational  grid 


Coordinates,  cell  number  Coordinates,  cell  number 


Figure  6.19:  Integral  radiation  heat  flux  depending  on 
order  of  used  quadrature,  for  10  spectral  groups  and 
NK=30 


Figure  6.20:  Integral  radiation  heat  flux  depending  on 
number  of  spectral  groups  for  S6  approximation  and 


NK=30 


6.5  Random  models  of  atomic  lines  [48] 

Radiation  heat  transfer  behind  shock  waves  in  air  is  considered.  The  problem  is  solved  in  one  dimensional 
(ID)  approach  for  temperature  profile  behind  the  shock  wave  fronts  which  is  shown  in  Figures  6.21,  a  at 
atmospheric  pressure.  Corresponding  distributions  of  the  molar  fractions  of  high  temperature  air  are 
shown  in  Figure  6.21,  d.  The  radiation  heat  transfer  equation  was  integrated  with  the  line-by-line  method 
and  with  the  random  model,  described  in  Section  4.1. 


For  the  analysis  two  spectral  ranges  Acd2,  Aco16  were  chosen  (see  Table  in  the  Figure  6.22).  The  range 
Acd2  contains  about  100  lines  of  atoms  N,  O,  formed  at  transitions  from  exited  energy  states.  In  the  each 
spectral  point  of  this  range  the  optical  thickness  is  appreciably  less  than  unit.  The  second  spectral  range 
Aco16  is  characterized  by  the  following  attributes:  the  spectral  range  contains  atomic  and  ionic  lines, 
formed  at  transitions  from  ground  and  lower  states;  the  optical  thickness  in  centers  of  the  strongest  lines 
surpasses  unit. 
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The  typical  distributions  of  integrated  radiation  fluxes  and  them  divergency  are  shown  in  Figures  6.21,b,c. 
Areas  of  positive  flux  divergency  correspond  primary  to  emission  of  radiation,  and  negative  ones 
corresponds  to  absorption.  The  observable  absorption  of  radiation  in  the  ultra  violet  part  of  the  spectrum  is 
caused  generally  by  the  photoionization  and  photodissociation,  and  also  by  some  electronic  bands  of 
diatomic  molecules. 


T ,  kK  FLUX,  Ns,tt/c=m+*2 


Figure  6.21 :  Distribution  of  temperature  and  total  radiation  characteristics  behind  shock  wave  with 
Tmax  =12000  K:  a)  temperature,  b)  one-side  (upper  and  lower  dashed  curves),  and  full  radiative  fluxes, 
c)  divergence  of  the  total  flux,  d)  species  molar  fraction 

Let  us  consider  radiative  characteristics  in  two  points  behind  the  shock  wave  front  with  coordinates 
xi  =0.5  and  x2  =3.7  cm.  In  the  first  of  the  two  points  the  spectral  and  integrated  fluxes  Mf  (xj)  give 
representation  about  radiation  of  the  layer  on  the  left  its  part.  Coordinate  x2  corresponds  to  right  boundary 
of  high-temperature  part  of  the  shock  wave.  The  values  of  spectral  and  total  fluxes  in  the  point  x\  are  of 
interest  for  the  purposes  of  the  analysis  of  photochemical  processes  in  the  cold  air.  The  values  of  radiation 
fluxes  in  the  point  x2  are  of  interest  for  the  analysis  of  processes  in  low-temperature  regions  of  the  shock 
wave. 


The  distributions  of  spectral  radiation  fluxes  in  the  points  x\  and  x2  are  shown  in  Figure  6.22.  On  can 
observe  some  spectral  regions,  in  which  the  leaving  radiation  is  close  to  radiation  of  the  absolutely  black 
body. 

Figures  6.23  and  6.24  show  spectral  absorption  coefficients  in  the  points  x\  and  x2,  and  also  spectral 
radiating  fluxes,  obtained  at  line-by-line  calculations  in  the  regions  Aco2  and  Aco16 .  It  is  seen  that  the 
atomic  lines  spectra  in  these  regions  differ  by  location,  intensities  and  half-widths. 

The  spectral  radiating  fluxes  Mf  (xx)  and  Mf  (x2)  in  the  range  Acd2  are  rather  close.  It  can  be  explained 
by  small  optical  thickness  in  all  its  spectral  points.  In  a  vicinity  of  the  strongest  lines  in  spectral  range 
Aco16  the  leaving  radiation  Mf  (xj)  almost  reaches  radiation  of  the  black  body.  Presence  of  the  same 
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lines  (low-broadening  in  low-temperature  part  of  the  plasma  layer)  results  in  strong  absorption  of  the 
spectral  radiation. 
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Figure  6.22:  Spectral  distribution  of  the  positive  radiation  flux  at  the  point  x=3.7cm,  and  spectral 
distribution  of  negative  radiation  flux  at  point  x=0.5cm;  solid  curves:  the  Planck  function. 
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Figure  6.23:  Spectral  absorption  coefficient  in  the  region  Aco2  =  10000-^13550  cm-1  at  the  point  x  =  3.7cm  (a) 
and  at  the  point  x  =  0.5cm  (c);  spectral  radiation  flux  in  the  region  Aco2  =10000h-13550  cm-1  at  the  point 
x  =  3.7cm  (b)  and  at  the  point  x  =  0.5cm(d).  7max=12000K 
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Figure  6.24:  Spectral  absorption  coefficient  in  the  region  Acg16  =  80400-^86800  cm-1  at  the  point  x  =  3.4cm  (a) 
and  at  the  point  x=0.5  cm  (c);  spectral  radiation  flux  in  the  region  Aco16  =80400-86800  cm-1  at  the  point 
x  =  3.4cm  (b)  and  at  the  point  x=0.5  cm  (d).  7"max  =12000  K 


To  give  general  representation  about  influence  of  atomic  and  ionic  lines  on  heat  radiative  transfer  in 
plasma  layers  the  distributions  of  total  radiation  fluxes  inside  the  layer  are  presented  in  Figure  6.25.  In 
spite  of  the  fact  that  in  separate  spectral  ranges  atomic  and  ionic  lines  give  increase  of  radiation  fluxes 
more  than  on  the  order,  total  radiation  grows  not  more  than  in  two-three  times.  Atomic  and  ionic  lines 
exert  largest  influence  on  the  radiative  fluxes  in  the  infra-red  area  of  the  spectrum,  where  the  optical 
thickness  of  the  layer  is  small.  In  the  ultra-violet  part  of  the  spectrum  the  influence  of  separate  groups  of 
strong  lines  is  not  very  great.  Either  the  share  of  these  sites  in  integrated  radiation  is  small,  or  the 
investigated  plasma  layers  here  have  optical  thickness  appreciably  larger  than  unit  and  integrated  fluxes 
poorly  differ  from  flows  of  an  absolutely  black  body. 
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Figure  6.25:  Distributions  of  the  one-side  (upper  and  lower  curves),  and  total  radiation  fluxes 
in  the  plane  layer  with  Tmax  =12000  K:  a)  atomic  lines  are  taken  into  account;  solid  curve  -  line-by-line 
calculation,  dashed  curve  -  group-random  model;  b)  without  atomic  lines 

Let  us  now  consider  results  of  estimations  of  radiation  heat  fluxes  obtained  with  the  use  of  group-random 
model.  At  preparation  of  group  functions  of  the  atomic  lines  parameters  the  following  principles  of  the 
grouping  were  used:  1)  inside  given  group  there  are  the  lines  only  one  species;  2)  the  lines  are  grouped  on 
sizes  of  constants  C4  (see  Section  4.2):  C4  =  0.001-0.01  for  the  first  group,  C4  =0.01-0.1  for  the  second 
group,  C4  =  0. 1-1.0  for  the  third  group,  C4  =  1.0-10  for  the  fourth  group. 

The  mathematical  expectation  of  the  oscillatory  strength  and  half-width  of  lines  inside  each  group  were 
calculated  as  the  average  arithmetic  of  their  sizes.  An  exponential  low  was  used  for  probability  densities 
of  oscillator  strengths,  and  the  delta-function  was  used  for  probability  densities  of  half-width  of  the  atomic 
lines. 

The  relative  discrepancy  between  the  line-by-line  and  random  numerical  simulation  results  are  calculated 
as  follows: 


M0  1  rm  ~  Mr)  i  7 _ 

e  =  X  + - ^100% , 

Moxis 

where  M^lls,MQlrm  are  the  radiation  fluxes,  obtained  at  the  line-by-line  integration  and  by  the  group- 

random  model.  Comparison  of  the  radiation  fluxes  with  and  without  atomic  lines,  as  well  as  the  radiation 
fluxes  obtained  by  the  line-by-line  and  the  random  model  are  shown  in  Figure  6.25.  It  is  possible  to 
conclude  that  the  use  of  the  group-random  model  allows  predict  total  radiation  fluxes  with  accuracy  not 
worse  than  ~30  %.  It  should  be  emphasized,  that  by  use  of  simple  statistical  model  (at  association  of  all 
lines  in  the  one  group  on  C4)  these  errors  can  reach  ~  200  %.  Comparison  of  the  radiation  fluxes  with  and 
without  atomic  lines  shows  significant  influence  of  the  atomic  lines  in  the  case  under  consideration. 

6.6  Macro-random  model  [20,  21] 

The  macro-random  model  of  vibrational  bands  of  molecules  H20  and  C02  was  applied  for  calculation  of 
radiative  characteristics  in  parallel-plane  inhomogeneous  layers.  Schematic  of  the  problem  of  radiation 
heat  transfer  is  shown  in  Figure  6.26.  Temperature  and  molar  fractions  of  gas  species  distributions  (C02, 
H20)  was  set  in  the  following  form: 
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T(x)  =  T0+(Tmax  -T0)exp 


x-xc 

n 

dx 

ck(x)  =  ck, 0  +  (cmax  -  ck,0  )exP 


X~Xn 


dv 


(57) 


where  7^ ,  T0 ,  ck  mdx ,  ck  0  are  the  maximal  and  minimal  values  of  temperature  and  concentration  of  the  k- 
th  species;  dx ,  n  are  the  parameters,  determining  a  degree  of  non-homogeneity  of  the  layer.  Figure  6.27 
shows  typical  distributions  of  temperature  and  molar  fractions  in  the  layer  T0  =300K,  Tmax  =2500K, 
(cH2O,0  =0.01,  CH2o,max  =0.7,  CCq2  0  =0.01,  Cc o2,max  =0.3). 


Figure  6.26:  Scheme  of  the  problem 

Calculations  of  radiation  heat  transfer  in  the  plane-parallel  layer  were  performed  with  the  “line-by-line” 
approach,  and  with  the  macro-random  model.  There  are  four  vibrational  bands  of  a  molecule  H20  and  four 
vibrational  bands  of  a  molecule  C02  were  taken  into  account  in  the  following  spectral  region 
AQ  =  1000-^10000  cm-1.  Parameters  of  vibrational  bands  were  borrowed  from  [22]. 

At  the  “line-by-line”  calculations  the  spectral  range  AQ  was  broken  on  elementary  spectral  regions  Acoz , 
in  the  centre  of  each  co,  the  total  absorption  coefficient  (averaged  on  rotational  structure)  was  calculated 
as  follows: 


gh2o  c 

k(©,.)  =  0.01  Z  77^exP 

g= 1  C3  ,g 


(Dj  -C0g 


C 


3  ,g 


S2  C, 


Ph2o  +0-01  Z  7T^exP 


g=l  C3,g 


C 


3,g 


Pco9 


(58) 


where  GH2q,  GCq2  are  the  numbers  of  vibrational  bands,  taken  into  account  for  H20  and  C02  accordingly; 
Ph20’Pco2  are  densities  of  H20  and  C02.  The  pressure  in  the  layer  was  supposed  as  atmospheric 
( p  =latm),  therefore  density  of  the  k- th  molecular  component  (with  molecular  weight  Mk )  is  calculated 
under  the  formula: 
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Pk  =MkxkP 


105 

8314 


kg/m3. 


The  absorption  coefficient  (58)  was  supposed  constant  in  the  limits  of  Aco^ . 


Spectral  absorption  coefficient  in  two  points  inside  the  plane-parallel  layer,  used  in  the  “line -by-line” 
calculations  is  shown  in  Figure  6.28.  Distributions  of  total  half-moment  functions  in  the  plane-parallel 
layer  are  shown  in  Figure  6.29.  Figure  6.29  shows  distribution  of  the  total  integrated  flux  obtained  in  the 
“line-by-line”  calculations: 


Wr(x)  =  M1+(x)  +  Mf  (jc), 

and  also  relative  error  of  these  calculation  with  the  use  of  the  macro-random  model  (WMSM  ( .v) ): 

K(*)-Wmw(x)| 

K(*)| 

Figure  6.30  shows  distribution  of  divergence  of  the  total  integrated  flux  and  error  of  its  calculation  with 
the  use  of  the  macro-random  model.  Divergence  of  radiating  flux  defines  volumetric  capacity  of  heat 
release,  therefore  this  function  is  first  of  all  necessary  for  the  solution  of  radiative-convective  interaction. 

One  can  see  that  the  error  of  prediction  of  the  radiation  function  with  the  macro-random  model  is  about 
~  30%,  that  is  quite  acceptable  for  the  practice  of  calculations  of  radiative  heat  transfer.  In  separate  very 
narrow  sites  this  error  can  reach  up  to  several  hundreds  of  percents.  However  it  is  observed  only  in  those 
points  inside  a  layer,  where  the  total  integrated  flux  or  divergence  of  the  flux  are  close  to  zero,  therefore 
the  specified  errors  do  not  worsen  the  accuracy  of  the  given  approximate  method. 


Figure  6.27:  Distribution  of  temperature  (solid  line,  scale  at  the  left),  and  molar  fractions  of  water  vapor  (long 
dashed  line;  scale  at  the  right)  and  carbon  dioxide  (short  dashed  line;  scale  at  the  right) 
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RESORPTION  COEFFICIENT,.  1/CM 


Figure  6.28:  Spectral  absorption  coefficient,  averaged  on  rotational  structure  in  points  of  the  layer  with 

coordinates  x  =  0  (1 )  and  x  =  H  (2) 


MOMENT  FUNCTIONS 


Figure  6.29:  Total  half-moment  characteristics  in  the  plane-parallel  layer 
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Figure  6.30:  Density  of  the  total  integrated  radiating  flux  in  the  layer  (solid  line;  scale  at  the  left)  and  error 
of  its  calculation  with  use  of  macro-random  model  (dashed  line;  scale  at  the  right) 


DIVERGENCE,  W/cm*+3  ERROR 


Figure  6.31 :  Divergence  of  the  total  integrated  flux  in  the  layer  (solid  line;  scale  at  the  left)  and  error 
of  its  calculation  with  the  use  of  macro-random  model  (dashed  line;  scale  at  the  right) 

6.7  The  Monte-Carlo  method 

6.7.1  Monte- Carlo  prediction  of  signature  of  model  solid  rocket  motor  plume  on  non-orthogonal 
grids  [31] 

Imitative  Monte-Carlo  algorithm  of  quasi-random  sampling  (see  Section  5.7)  is  used  here  for  prediction  of 
spectral  signature  of  multi-phase  plume  of  model  solid  rocket  motor  (SRM).  Gas  dynamic  functions  of  the 
model  SRM  plume  are  shown  in  Figures  6.32.  These  are:  gas  temperature  ( a ),  temperature  of  solid 
particles  (/?),  volume  concentration  of  the  monodisperse  cloud  of  A1203  particles  with  radius  of  particles  of 
2  microns  (c),  molar  fractions  of  H20  ( d)  and  C02  ( e ).  Non-orthogonal  computational  grid  used  for  these 
calculations  is  shown  in  Figure  6.32,  f. 
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X,  cm 


X,  cm  X,  cm 


X,  cm  X,  cm 


R,  cm  R,  cm  R,  cm  R> cm  R,  cm  R,  cm 


Figure  6.32:  Gas  temperature  (a),  temperature  of  AI2O3  particles  ( b )  and  volume  concentration  of  AI2O3 
particles  in  model  scattering  plume  of  SRM  (c)  on  the  non-orthogonal  axysimmetric  grid 


Calculations  were  performed  for  infrared  region  of  spectrum  corresponding  to  center  of  strong  absorption 
band  of  water  vapors  (the  2.7p-band),  Aco  =  3600  -r-  3620  cm1.  Pressure  in  the  plume  was  assumed  to  be 
p  =  1  atm.  A  spectral  absorption  coefficient  in  spectrum  of  the  Lorentzian  rotational  lines  inside  and 
outside  of  not  scattering  and  scattering  model  SRM  plumes  is  shown  in  Figure  6.33. 


Wavenumber,  1 /cm  ^  Wavenumber,  1 /cm  ^ 

Figure  6.33:  Spectral  absorption  coefficient  in  two  points  inside  scattering  model  plume;  a:  solid  line  -  X=  0, 
R  =  0,  dashed  line  -  X=  0,  R=  10  cm;  b :  solid  line  -  X=  600  cm,  R=  0,  dashed  line  -  X=  600  cm,  R=  30  cm. 
Albedo  of  the  single  scattering  (the  bold  solid  line)  inside  scattering  model  plume  (a:  X=  R-  0;  b:  X=  600  cm, 

R-  30  cm) 
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The  computational  case  of  the  light-scattering  plume  is  characterized  by  presence  of  spectral  regions  with 
very  high  single-scattering  albedo,  that  is  by  the  very  high  probability  of  scattering  (Figures  6.33,a,b). 

Results  of  the  Monte-Carlo  calculations  of  integral  and  spectral  signatures  of  model  non-scattering  plume 
on  the  non-orthogonal  grid  are  shown  in  Figures  6.34-6.36.  Dependences  of  integral  directional  signature 
of  the  plume  from  the  angle  of  observation,  calculated  with  use  of  three  described  imitative  methods,  are 
shown  in  these  figures.  The  Monte-Carlo  simulation  algorithms,  which  use  methods  of  the  quasi-random 
sampling  (Figure  6.34)  and  the  Maximum  Cross-Section  (Figure  6.35),  give  close  results. 


Spectral  directional  emissivity  (Spectral  signature),  W/(ster*mcm) 


Integral  directional  emissivity  (Integral),  W/ster 


Figure  6.34:  Left:  Integral  signature  of  non-scattering  plume  vs.  angle  0  of  the  plume  observation.  Right: 
Spectral  signature  of  SRM  in  two  directions  to  X-axis  (0  =  90°  and  0  =  150°).  The  direct  LBL  Monte  Carlo 
method  on  non-orthogonal  axysimmetric  grid  (number  of  points  along  photon  trajectory:  30) 


Integral  directional  emissivity  (Integral),  W/ster 


Angle,  degree 


Spectral  directional  emissivity  (Spectral  signature),  W/(ster*mcm) 


Wavenumber,  1/cm 


Figure  6.35:  Left:  Integral  signature  of  non-scattering  model  plume  of  SRM  vs.  angle  0  of  the  plume 
observation.  Right:  Spectral  signature  of  non-scattering  model  plume  of  SRM  in  two  directions  of  observation 
to  X-axis  (0  =  90°  and  0  =  150°).  The  LBL  Monte  Carlo  method  with  the  Maximal  Cross-Section  algorithm  on 

non-orthogonal  axysimmetric  grid 
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The  hybrid  method  (Figure  6.36)  gives  the  overestimated  results  in  the  normal  direction  to  X-axis. 
However,  it  is  necessary  to  note,  that  the  optimization  of  approximating  parameters  of  molecular  rotational 
lines  statistical  models  was  not  performed  in  this  case. 

Integral  directional  emissivity  (Integral),  W/ster  Integral  directional  emissivity  (Integral),  W/ster 


Figure  6.36:  Integral  signature  of  non-scattering  model  plume  of  SRM  vs.  angle  0  of  the  plume  observation 
(the  angle  6  is  counted  off  from  the  X-axis).  Left:  the  Hybrid  Monte  Carlo  method.  Right:  the  Direct  LBL  Monte 
Carlo  method  (number  of  points  along  photon  trajectory:  30) 

The  spectral  dependences  of  signatures  at  two  angles  of  observation,  which  are  presented  in  Figures  6.34 
and  6.35,  show  some  differences  in  separate  narrow  spectral  ranges.  These  differences  can  be  referred  to 
statistical  errors  of  the  Monte-Carlo  algorithms. 


Calculation  results  for  integral  and  spectral  signatures  for  scattering  plume  of  the  model  SRM  are  shown 
in  Figures  6.37-6.39. 


Integral  directional  emissivity  (Integral),  W/ster 


Spectral  directional  emissivity  (Spectral  signature),  W/(ster*mcm) 
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Fig.6.37  Integral  signature  (left)  of  scattering  model  plume  of  SRM  vs.  angle  0  of  the  plume  observation. 
Spectral  signature  (b)  of  scattering  model  plume  of  SRM  in  two  directions  to  X-axis  (0  =  90°  and  0  =  150°).  The 
direct  LBL  Monte  Carlo  method  (number  of  points  along  photon  trajectory:  30) 
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As  well  as  for  the  non-scattering  medium,  the  results  obtained  by  two  LBL-methods  are  rather  close.  It 


concerns  both  to  the  total  signature  (Figures  6.37, 
6.37,b  and  6.38,b). 


Integral  directional  emissivity  (Integral),  W/ster 


and  6.37,b)  and  to  the  spectral  signature  (Figures 

Spectral  directional  emissivity  (Spectral  signature),  W/(ster*mcm) 


Wavenumber,  1/cm 


Figure  6.38:  Integral  signature  (a)  of  scattering  model  plume  of  SRM  vs.  angle  0  of  the  plume  observation. 
Spectral  signature  of  scattering  model  plume  of  SRM  in  two  directions  to  X-axis  (0  =  90°  and  0  =  150°).  The 
LBL  Monte  Carlo  method  with  the  Maximal  Cross-Section  algorithm 


The  hybrid  method  also  shows  rather  good  accuracy  in  the  case  under  consideration  (Figure  6.39). 


Integral  directional  emissivity  (Integral),  W/ster 


Figure  6.39:  Integral  signature  of  scattering  model  plume  of  SRM  vs.  angle  0  of  the  plume  observation. 

The  Hybrid  Monte  Carlo  method 

Presented  data  allow  make  conclusion  about  effect  of  scattering  on  a  signature  of  plumes.  Let  us  compare, 
for  example,  the  spectral  signature  calculated  by  the  LBL-method  with  the  Maximum  Cross-Section 
method  for  non-scattering  (Figure  6.35)  and  for  scattering  (Figure  6.38)  plumes.  Scattering  processes 
become  apparent  appreciably  for  angles  of  observation  far  from  90.  This  is  because  in  paraxial  directions 
the  absorption  probability  of  photons  increases.  Integral  radiation  in  the  normal  direction  to  the  axis  of 
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symmetry  in  this  case  changes  insignificantly,  however  the  spectral  distribution  of  radiation  undergoes 
noticeable  changes:  the  spectral  signature  of  the  light-scattering  plume  has  smaller  oscillations  on  the 
spectrum.  This  fact  also  can  be  explained  by  the  presence  of  additional  absorber  (particles  of  A1203)  and 
by  much  higher  probability  of  photon  survival  at  any  elementary  collision  act  of  photons  with  particles  of 
the  medium. 

6.7.2  Monte-Carlo  prediction  of  spectral  emissivity  of  three-dimensional  plumes  [33] 

Developed  algorithm  of  imitating  modeling  of  radiation  heat  transfer  in  a  three-dimensional  case  was  used 
for  analysis  of  spectral  and  integral  signatures  of  model  jet  shown  in  Figures  6.40-6.42.  Two  calculation 
cases  were  analyzed.  The  first  one  is  the  weak  light-scattering  plume  (concentration  of  A1203  particles  in 
the  exit  cross  section  of  each  SRM  is  np  =  5  x  103  cm'3),  and  the  second  is  the  strong  light-scattering  plume 
(concentration  of  A1203  particles  A1203  exit  cross  section  of  each  SRM  is  np  =  107cnT3).  The  average 
radius  of  particles  was  presumed  equal  to  2  and  20  microns. 

Distribution  of  temperatures  of  gas  and  condensed  phase  were  set  on  the  basis  of  the  approached  analytical 
solution  of  the  problem  about  two-phase  plume,  and  distribution  of  mass  fractions  of  optically  active 
components  was  set  on  the  base  on  theory  of  analogy  of  heat  and  mass  transfer  processes.  Spectral  optical 
properties  of  combustion  products  were  calculated  with  use  of  the  computer  system  ASTEROID  [1,  47]. 
Calculations  of  the  optical  properties  of  particles  A1203  were  performed  with  the  use  of  the  Mie  theory 
[29].  Spectral  calculations  of  radiation  heat  transfer  were  performed  for  100-group  model  of  optical 
properties  with  averaging  of  rotational  structure  of  molecular  spectrum.  Trajectories  of  105  groups  of 
photons  were  modeled  in  each  spectral  group.  Control  calculations  were  executed  at  modeling  of 
106  groups  of  photons  in  each  spectral  group. 
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Figure  6.40:  Gas  temperature  in  three-block  plume 


RTO-EN-AVT-1 62 


10-®$ 


Radiation  Modeling  in  Shock-Tubes  and  Entry  Flows 


\ 


TS 

1.92E+03 

1.85E+03 

1.77E+03 

1.69E+03 

1.62E+03 

1.54E+03 

1.46E+03 

1.38E+03 

1.31  E+03 
1.23E+03 
1.15E+03 
1.08E+03 
1.00E+03 
9.23E+02 
8.46E+02 
7.69E+02 
6.92E+02 
6.15E+02 
5.38E+02 
4.62E+02 
3.85E+02 
3.08E+02 

2.31  E+02 
1.54E+02 
7.69E+01 


! 

I 


■ 


Ysl 

8.55E+06 
5.33E+06 
3.32E+06 
2.07E+06 
1.29E+06 
8.03E+05 
5.00E+05 
3.12E+05 
1.94E+05 
1.21E+05 
7.53E+04 
4.69E+04 
2.92E+04 
1.82E+04 
1.14E+04 
7.07E+03 
4.41  E+03 
2.74E+03 
1.71  E+03 
1.07E+03 
6.64E+02 
4.14E+02 
2.58E+02 
1.61  E+02 
1.00E+02 


Figure  6.41 :  Temperature  of  AI2O3  particles  in  three- 
block  plume 


Figure  6.42:  Concentration  of  AI2O3  particles  in  three- 
block  plume 
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Figure  6.43:  Mass  fraction  of  H20 


Figure  6.44:  Mass  fraction  of  CO2 
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Signatures  of  radiation  of  three -block  plume  are  shown  in  Figures  6.45  and  6.46.  Created  calculation 
algorithm  demonstrates  symmetry  of  received  results  in  relation  to  turn  on  90°.  In  the  first  case  (Figure 
5.45)  signature  of  the  plume  in  YZ  plane  exceeds  radiating  ability  from  XZ  plane.  A  reverse  situation  is 
observed  in  the  second  case  (Figure  5.46).  Marked  regularities  for  integral  radiation  are  repeated  also  for 
the  spectral  signature.  In  the  considered  case,  as  before,  radiating  ability  is  formed  in  the  centers  of 
vibrational  bands,  and  the  solid  phase  radiation  is  distinctly  visible  in  transparency  windows  of  vibrational 
bands. 


Integral  directional  emissivity  (Integral),  W/ster  Integral  directional  emissivity  (Integral),  W/ster 


Figure  6.45:  Integral  signature  of  the  three-block  plume  with  high  level  of  scattering; 

Observation  of  plane  at  x=xmax  (left)  and  plane  y  =  ymax  (right) 

Spectral  directional  emissivity  (Spectral  signature),  W/(ster*mcm)  Spectral  directional  emissivity  (Spectral  signature),  W/(ster*mcm) 


Wavenumber,  1  /cm  Wavenumber,  1  /cm 

Figure  6.46:  Spectral  signature  of  the  three-block  plume  with  high  level  of  scattering;  Observation  of  plane  at 

*=*max  (left)  and  plane  y  =  ymax  (right) 
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6.7.3  Monte- Carlo  simulation  with  line-by-line  spectral  modelling  [49] 

Computing  code  Plume-3D-MC  was  used  for  interpretation  of  the  ERIM  experimental  data  [2],  and  to 
systematic  study  of  selective  radiation  heat  transfer  from  cloud  of  hot  gases  through  cold  air  atmosphere 
[50].  To  focus  attention  on  verification  of  spectral  optical  models  included  into  code  Plume-3D-MC  the 
problem  of  prediction  of  spatial  radiation  characteristics  was  significantly  simplified  in  comparison  with 
previous  study.  In  this  case  homogeneous  cylindrical  volume  of  height  Hhc  =  60  cm  and  radius 
Rhc  =  30  cm  contains  mixture  of  water  vapour  or  carbon  dioxide  with  molecular  nitrogen  at  temperature 
T  ~  1200  K  and  pressure  p  =  0.1  atm.  Such  conditions  were  experimentally  studied  in  [50]. 


Cylindrical  volume  involving  heated  mixture  of  H20/N2  at  temperature  T  =  1202K  and  total  pressure 
p  =  0.1  was  considered  in  the  second  case.  Molar  fractions  of  gas  species  are  xHlo  =0.5,  x^2  =0.5. 

Numerical  prediction  of  spectral  radiation  intensity  of  the  heated  volume  in  normal  direction  to  plane 
surface  are  shown  in  Figure  6.47.  Optical  properties  of  heated  H20  were  borrowed  from  NASA  standard 
infrared  radiation  model  [35],  where  rotational  line  structure  was  averaged  in  each  spectral  group  of 
Aco  =  25cm“1.  Experimental  data  [50]  are  presented  for  spectral  region  AQ  =  3000-^4200  cm'1.  These 
calculations  were  performed  with  the  same  spectral  resolution.  Spectral  absorption  coefficient  calculated 
at  these  conditions  is  shown  in  Figure  6.48. 


Group  directional  emissivity  (Group  signature),  W/fste^cm'^cm2) 


3000  3200  3400  3600  3800  4000  4200 


Wavenumber,  1/cm 


Wavenumber,  1/cm 


FIGURE  6.47:  TEST  9R.  ERIM  EXPERIMENTAL  DATA  FIGURE  6.48:  SPECTRAL  ABSORPTION  COEFFICIENT 
FOR  H20  HOT  CELL  RADIANCE  (UPPER  DOTTED  OF  HEATED  H20  FOR  CONDITIONS  OF  TEST  9R. 

LINE),  HOT-THROUGH-COLD  RADIANCE  (LOWER  NASA  STANDARD  INFRARED  RADIATION  MODEL 

DOTTED  LINE),  AND  PLUME-3D-MC  NUMERICAL 
PREDICTION  (SOLID  LINE)  OF  THE  HOT  CELL 
RADIANCE.  NASA  STANDARD  INFRARED  RADIATION 
MODEL;  AVERAGED  ROTATIONAL  LINE  STRUCTURE; 

SPECTRAL  REGION  FOR  EACH  SPECTRAL  GROUP  IS 
A(o  =  25  CM-1 


Observed  in  Figure  6.47  good  agreement  between  theoretical  prediction  and  ERIM  experimental  data 
indicates  that  optical  model  [35]  provides  adequate  data  for  the  case  under  consideration. 
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Cylindrical  volume  involving  heated  mixture  of  C02/N2  at  temperature  T  =  1200K  and  total  pressure 
p  =  0.1  was  considered  in  the  second  case.  Molar  fractions  of  the  gas  species  are  Xq2 0  =0.5,  xNl  =0.5  . 
Other  initial  conditions  were  the  same  as  in  the  first  calculation  series. 


Numerical  simulation  results  for  spectral  intensity  obtained  with  the  use  of  code  Plume-3D-MC  are  shown 
in  Figure  6.49.  The  NASA  standard  infrared  radiation  model  [35]  was  used  in  these  calculations. 
Averaged  over  rotational  line  structure  spectral  absorption  coefficient  of  hot  mixture  C02/N2  for 
conditions  of  experimental  research  [50]  (Test  5)  is  shown  in  Figure  6.50. 

Considered  above  numerical  prediction  of  spectral  radiation  intensity  of  heated  volumes  of  H20/N2  and 
C02/N2  for  conditions  of  experiments  [50]  (Test  9R  and  Test  5),  were  obtained  with  the  data  accumulated 
in  NASA  standard  infrared  model  [35].  It  is  of  great  practical  interest  to  consider  possibility  to  use  other 
optical  models  for  interpretation  of  the  experimental  data.  Results  of  such  investigations  are  presented  in 
[49]. 


2150  2200  2250  2300  2350  2400  2450 

Wavenumber,  1/cm 


Wavenumber,  1/cm 


FIGURE  6.49:  TEST  5.  ERIM  EXPERIMENTAL  DATA  FIGURE  6.50:  SPECTRAL  ABSORPTION 

FOR  C02  HOT  CELL  RADIANCE  (DOTTED  LINE)  AND  COEFFICIENT  OF  HEATED  C02  AT  CONDITIONS  OF 

PLUME-2D-MC  NUMERICAL  PREDICTION  (SOLID  TEST  5;  NASA  STANDARD  INFRARED  RADIATION 
LINE)  OF  THE  HOT  CELL  RADIANCE.  THE  NASA  MODEL  [35] 

STANDARD  INFRARED  RADIATION  MODEL; 

AVERAGED  ROTATIONAL  LINE  STRUCTURE; 

SPECTRAL  REGION  FOR  EACH  SPECTRAL  GROUP 

IS  Aro  =  20CM-1 


Next  calculation  case  corresponded  to  experimental  data  on  observation  of  hot  volumes  of  H20  from  large 
distances.  Schematic  of  the  problem  is  shown  in  Figure  6.51.  Initial  conditions  were  identical  to  real 
experimental  investigation.  Parameters  of  gas  mixture  inside  heated  volume:  p  =  0.1  is  the  total  pressure, 

atm;  /?h2o  =  0.05 ,  p^2  =  0.05  are  the  partial  pressures  of  species,  atm;  H  =  60  is  the  height  of  the  heated 
cylindrical  volume,  cm;  T  =  1202  is  the  temperature  of  gas  mixture,  K.  Parameters  of  surrounding  cold 
gas  are:  p  =  0.07  is  the  total  pressure,  atm;  Ph2o  =0.001,  p^2  =0.069  are  the  partial  pressures  of 
species,  atm;  L  =  10000  is  the  length  of  the  optical  path  of  observation,  cm;  T  =  293  is  the  temperature  in 
surrounding  gas,  K.  Note,  that  parameters  of  hot  gas  model  conditions  in  exhaust  rocket  plumes  at  altitude 
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about  ~16  km,  and  parameters  of  cold  gas  model  optical  path  of  ~  28.8  km  at  altitude  18  km.  Note  that  the 
calculation  domain  contains  to  sub-domains  (Figure  6.51).  The  first  one  exactly  corresponds  to  conditions 
of  experiment  Test-9R.  The  second  one,  which  models  volume  of  cold  gas,  is  100  times  smaller  than  in 
real  experimental  situation.  To  model  radiation  transfer  in  real  conditions,  the  following  approach  was 
used  in  code  Plume-3D-MC.  An  effective  pressure  of  H20  in  surrounding  gas  was  introduced: 

pk2o  =  100Ph2o  =  0.1  atm,  that  is  xHlQ  =  Ph2o  / P  =  1-429  .  But  this  effective  pressure  can  be  used  only  for 
prediction  of  averaged  absorption  coefficient  (parameter  S/d)  for  narrow-band  random  model  of 
rotational  lines.  To  calculate  broadening  of  rotational  lines  in  real  surrounding  conditions  (parameter 
y/d)  one  should  use  the  real  pressure  p  =  0.07  atm. 


Z 


Figure  6.51:  Schematic  of  the  3D  “Hot/Cold”  problem  of  visibility  of  heated  cylindrical  volume 
of  H20/N2  and  CO2/N2  through  cold  atmosphere 

Prediction  of  hot  cell  and  hot-thought-cold  spectral  radiance  by  code  Plume-3D-MC  and  NASA  standard 
infra-red  model  with  averaged  rotational  line  structure  is  shown  in  Figure  6.52  and  6.53.  One  can  see  from 
these  data  that  essentially,  the  problem  is  as  follows:  a  major  portion  of  heat  radiation  is  absorbed  in 
surrounding  gas,  containing  resonant  absorbing  cold  gas.  Figure  6.51  shows  comparison  of  experimental 
data  and  numerical  prediction  of  hot-thought-cold  spectral  radiance  for  conditions  of  Test-9R  experiments. 
Reference  to  Figure  6.51  shows  that  maximal  discrepancy  between  experimental  data  and  numerical 
prediction  is  observed  within  the  region  of  maximal  absorption  ability  of  H20  vibrational  bands.  At  the 
same  time,  it  is  well  known,  that  real  spectra  (vibrational  bands)  of  emission  and  absorption  contain  large 
number  of  molecular  rotational  lines,  therefore  it  might  be  assumed  that  the  reason  for  the  large 
discrepancy  between  experimental  and  numerical  data  is  in  using  of  non-adequate  spectral  optical  model. 
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Figures  6.54  and  6.55  show  numerical  prediction  of  hot-through-cold  spectral  radiance  for  ERIM 
experimental  conditions  Test-9R  (mixture  H20/N2),  obtained  by  code  Plume-3-MC  with  the  use  of 
narrow-random  models  [22].  Numerical  simulation  data  presented  in  Figures  6.54  and  6.55  correspond  to 
NASA  standard  infrared  radiation  optical  model  [35].  The  following  two  spectral  optical  models  of 
radiation  transfer  were  used  for  both  databases: 

•  the  narrow-band  random  model  for  optically  thin  layers  (the  method  of  smoothing  coefficients; 
see  Section  5.3); 

•  the  narrow-band  random  model  of  real  spectrum  (Figure  6.55). 


Wavenumber,  1/cm 


3200  3400  3600  3800  4000  4200 


Wavenumber,  1/cm 


FIGURE  6.52:  SPECTRAL  SIGNATURE  OF 
HOMOGENEOUS  CYLINDRICAL  VOLUME  OF  HEATED 
H20  UP  TO  T  =  1202  K; 

3D  CALCULATION  DOMAIN;  Nph  =104  ;  THE  MCLEDE 
ALGORITHM.  SOLID  LINE-  PREDICTION  OF  HOT  CELL 
RADIANCE  (TEST  9R);  DOTTED  LINE  -  PREDICTION  OF 
HOT-THROUGH-COLD  CELL  RADIANCE  (TEST-9R).  THE 
NASA  STANDARD  INFRARED  MODEL 


FIGURE  6.53:  SPECTRAL  SIGNATURE  OF 
HOMOGENEOUS  CYLINDRICAL  VOLUME  OF 
HEATED  H20  UP  TO  7  =  1202  K;  3D  CALCULATION 
DOMAIN;  Nph  =104  ;  THE  MCLEDE  ALGORITHM. 
CODE  PLUME-3D-MC.  SOLID  LINE  -  PREDICTION 
OF  HOT-THROUGH-COLD  CELL  RADIANCE, 
DOTTED  LINE  -  ERIM  EXPERIMENTS  FOR  TEST- 
9R 
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3000  3200  3400  3600  3800  4000  4200 

Wavenumber,  1/cm 

FIGURE  6.54:  TEST  9R.  ERIM  EXPERIMENTAL  DATA 
FOR  H20  HOT  CELL  RADIANCE  (UPPER  DOTTED 
LINE),  HOT-THROUGH-COLD  RADIANCE  (LOWER 
DOTTED  LINE),  AND  PLUME-3D-MC  NUMERICAL 
PREDICTION  (SOLID  LINE)  OF  THE  HOT  -THROUGH- 
COLD  RADIANCE.  THE  NASA  STANDARD  INFRARED 
RADIATION  MODEL;  THE  METHOD  OF  SMOOTHING 
COEFFICIENT;  SPECTRAL  REGION  FOR  EACH 
SPECTRAL  GROUP  IS 
Aco  =  25  CM-1 


3000  3200  3400  3600  3800  4000  4200 


Wavenumber,  1/cm 


FIGURE  6.55:  TEST  9R.  ERIM  EXPERIMENTAL  DATA 
FOR  H20  HOT  CELL  RADIANCE  (UPPER  DOTTED 
LINE),  HOT-THROUGH-COLD  RADIANCE  (LOWER 
DOTTED  LINE),  AND  PLUME-2D-MC  NUMERICAL 
PREDICTION  (SOLID  LINE)  OF  THE  HOT  -THROUGH- 
COLD  RADIANCE.  THE  NASA  STANDARD  INFRARED 
RADIATION  MODEL;  RANDOM  MODEL  (JLBL  =  2); 
SPECTRAL  REGION  FOR  EACH  SPECTRAL  GROUP 
IS  Aoo  =  25  CM-1 


With  reference  to  presented  numerical  simulation  data  (Figure  6.55),  it  can  be  seen  that  narrow-band 
random  models  provide  reasonable  prediction  of  hot-through-cold  spectral  radiance  for  ERIM 
experimental  conditions.  However,  some  discrepancies  between  predicted  and  experimental  data  (inside 
separate  spectral  regions)  have  engaged  our  attention,  and  it  substantiates  necessity  in  further  development 
of  spectral  optical  models  of  heated  gases  and  radiation  heat  transfer. 

And  finally,  prediction  of  hot-through-cold  radiance  with  the  use  of  line-by-line  spectral  models  based  on 
HITRAN-like  databases  will  be  considered.  Figures  6.56,a-d  show  numerical  prediction  of  hot  cell 
spectral  radiance  for  ERIM  experimental  conditions  Test-9R  (mixture  H20/N2),  obtained  by  code  Plume  - 
3D-MC  with  the  use  of  line-by-line  model.  These  numerical  simulation  data  were  obtained  with  the  use  of 
parameters  of  rotational  lines  from  HITRAN  database  [36]. 
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Wavenumber,  1/cm 


a 


b 


Wavenumber,  1/cm  Wavenumber,  1/cm 


C 


d 


Figure  6.56:  Line-by-line  prediction  of  spectral  intensity  for  conditions  of  ERIM  experiments  (Test  9R); 
Hot  cell  H20/N2  radiance;  spectral  resolution  is  Aco  =  0.0083  cm"1 


Figure  6. 56, a  shows  spectral  intensity  of  hot  cell  radiance  in  spectral  region  AQ  =  3000-^4400  cm'1.  This 
spectral  region  was  divided  by  56  spectral  groups.  The  Monte-Carlo  simulation  of  radiation  transfer  was 
performed  in  3000  spectral  points  inside  each  spectral  group.  Thus,  the  spectral  resolution  of  these 
calculations  was  about  ~  Aco  =  0.0083  cm'1.  Figures  6.56,b-d  show  the  spectral  intensity  inside 
successfully  decreasing  spectral  regions  (up  to  AQ  =  3450  -^3452  cm'1).  Figures  6.57,a-c  show  numerical 
prediction  of  hot-through-cold  spectral  radiance  for  ERIM  experimental  conditions  Test-9R  (mixture 
H20/N2),  obtained  by  code  Plume-3D-MC  with  the  use  of  line-by-line  model.  These  numerical  simulation 
data  were  also  obtained  with  the  use  of  parameters  of  rotational  lines  from  HITRAN  database  [36]. 


RTO-EN-AVT-1 62 


10-911 


Radiation  Modeling  in  Shock-Tubes  and  Entry  Flows 


0.0005 


0.0004 


0.0003 


0.0002 


0.0001 


- 

Nphot/Giuup=1000000 

- 

- 

From  upper  surface 

- 

WWW  ■ 

—  O--  ERIM  data  Test  9R  (Hot) 

ERIM  data  Test  9R  (Hot/Cold) 

- 

- 

- 

- 

IIP  Pli  WI 

% 

3200  3400  3600  3800  4000  4200 

Wavenumber,  1/cm 


a 


Group  directional  emissivity  (Group  signature),  W/(ster*cm  ^cm2) 


Wavenumber,  1/cm 

b 


Group  directional  emissivity  (Group  signature),  W/(ster*cm  ^cm2) 


C 

Figure  6.57:  Line-by-line  prediction  of  spectral  intensity  for  conditions  of  ERIM  experiments  (Test  9R); 
Hot-through-cold  radiance  (H20/N2);  spectral  resolution  is  Aco= 0.0083  cm'1 

With  reference  to  Figures  6.56  and  6.57  it  can  be  seen  that  the  spectral  intensity  is  changed  inside  range  of 
~5  orders.  Therefore,  it  is  easy  to  understand  that  the  problem  of  averaging  of  the  intensity  in  spectral 
region  of  ~  Aoig  ~  25  cm"1  is  very  complex  and  difficult  problem,  which  is  very  sensitive  to  many  input 
factors  of  used  models. 

6.7.4  Spectral  signatures  of  fireballs  generated  at  chemical  explosions  [51] 

Radiative  gas  dynamics  and  emission  of  heat  radiation  of  fireballs  generated  at  large-scale  accident  at 
chemical  industry,  natural  gas  transportation  or  rocketry  is  considered  in  the  section.  The  problem  will  be 
illustrated  on  the  example  of  rocket  accident  at  active  part  of  trajectory.  It  was  shown  at  qualitative 
analysis  of  processes  accompanying  explosion  of  rockets  on  a  launching  pad  or  in  flight  [52]  that  after  end 
of  shock-wave  stage  of  the  process,  formation  of  fireball  is  rather  probable.  The  fireball  is  a  burning  cloud 
of  a  mix  of  components  of  rocket  fuel,  which  emerges  in  the  environment  atmosphere  under  action  of  the 
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Archimedean  force,  involving  in  the  movement  big  masses  of  ambient  air.  The  process  of  the  fireball 
formation,  in  turn,  can  be  presented  as  the  following  two-stage  process: 

1)  Expansion  of  the  burning  mix  during  the  first  3+7  seconds  of  the  process .  Feature  of  this  phase  of 
the  process  is  the  fast  increasing  of  the  fireball  sizes.  Average  speed  of  such  expansion  reaches  to 
20  m/s.  Thermal  radiation  of  the  fireball  accrues  most  quickly  at  this  stage,  reaching  the  maximum 
at  5+1  second  of  the  process.  Then  intensity  of  heat  radiation  gradually  decreases  owing  to  intensive 
cooling  of  the  fireball. 

2)  Emersion  of  the  fireball  in  the  environment  atmosphere.  Feature  of  this  stage  is  fast  emersion  of  the 
fireball,  accompanying  with  involving  in  movement  of  the  big  masses  of  air  and  formation  of  large- 
scale  vortical  movement  in  the  free  atmosphere.  The  fireball  intensively  exchanges  heat  with 
ambient  air  during  its  emersion.  In  spite  of  the  fact,  that  its  thermal  radiation  is  not  so  high  (due  to 
relative  low  temperatures,  ~  500-^1000  K),  heat  exchange  by  radiation  plays  an  important  role  in 
the  fireball  dynamics. 

Multi-component  model  of  chemically  reacting  gas  generated  at  burning  of  rocket  fuel  (H20,  C02,  CO, 
H2)  was  suggested  and  studied  in  [51]  for  prediction  of  the  fireballs  dynamics  in  view  of  heat  exchange  by 
heat  radiation  and  its  influence  on  gas-dynamic  processes  (in  this  sense  the  term  radiating  gas  dynamics 
was  introduced  in  [52]).  It  will  be  shown  here  how  the  spectral  signature  registered  from  large  distance  of 
the  fireball  reflects  physical-chemical  processes  accompanying  the  phenomenon  under  consideration. 

Figures  6.58,  6.59  show  temperature,  mass  fraction  of  H20  and  C02  molecules,  and  velocity  vector  field 
in  the  fireball  at  consecutive  time  moments  after  a  rocket  explosion. 

Being  based  on  these  data  we  shall  describe  a  qualitative  picture  of  a  process  of  the  fireball  dynamics: 

1)  Rather  fast  expansion  of  the  fireball  (~3  s)  with  its  simultaneous  cooling  approximately  up  to 
~  1000  K  is  observed.  The  fireball  form  remains  practically  spherical  on  this  time  interval. 

2)  Vortical  movement  of  gas  is  formed  inside  and  outside  of  the  fireball  on  the  5th  second  of  the 
process. 

3)  To  the  same  time  moment  the  fireball  begins  to  lose  its  spherical  form  (Figure  6.58, a).  It  concerns 
both  to  temperature  distribution,  and  to  distribution  of  mass  fractions  of  molecular  species.  The  top 
boundary  of  the  fireball  at  this  time  moment  is  displaced  approximately  on  200  m,  reaching  height 
of  800  m. 

4)  Further,  the  fireball  continues  to  move,  being  increased  in  the  sizes  that  occur  because  of  involving 
in  the  movement  of  atmospheric  air,  which  is  heated  and  mixed  with  products  of  combustion  of  the 
rocket  fuel.  For  example,  through  1 1  cex  after  explosion,  the  top  boundary  of  the  fireball  heated  up 
area  reaches  height  of  1000  m  (Figure  6.58,b),  and  to  21  s  -  1200  m  (Figure  6.58,d).  From  these 
figures  the  increase  of  the  fireball  radial  sizes  is  clearly  visible. 

5)  As  it  was  already  marked,  the  cooling  of  the  fireball  is  observed  in  process  of  its  rise.  A  rate  of  the 
cooling  is  the  most  quick  at  the  first  5  s,  due  to  significant  radiating  losses  of  energy.  Gradually  the 
basic  mechanism  of  the  energy  loss  due  to  thermal  radiation  is  replaced  by  the  convective  heat 
exchange,  and  the  rate  of  the  cooling  is  reduced.  By  the  21st  second  of  the  process  the  temperature 
inside  the  fireball  falls  approximately  to  600  K. 

6)  It  is  visible  from  Figure  6.58  that  to  21st  second  of  the  process  the  vortical  movement  of 
atmospheric  air  and  burning  products  covers  region  more  than  1  km  on  height,  and  ~  500  -f  600  m  in 
radial  direction. 
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Figure  6.58:  Temperature  distribution  in  fireball  at  consecutive  time  moments: 
(a)  t  =  5  s,  (b)  t  =  11  s,  (c)  t  =  15  s,  (d)  t  =  21  s 
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Figure  6.59:  Velocity  field  in  fireball  at  consecutive  time  moments:  (a)  t  =  5  s,  (b)  t  =  11  s,  (c)  t  =  15  s, 

(d)  t  =  21  s;  velocity  scale  (50  m/s)  is  shown  in  right  top  corner 

Strong  influence  of  radiating  processes  on  dynamics  of  fireballs  was  established  in  [52].  It  means  that  in 
spite  on  the  fact  that  temperature  inside  fireball  is  not  very  high,  the  radiation  emission  and  reabsorption 
have  significant  influence  on  the  fireball  dynamics  due  to  its  large  sizes.  For  example,  the  full  radiative 
gas  dynamic  model  gives  maximum  temperature  inside  the  fireball  ~  1140  K  to  the  5th  second  after 
explosion.  If  we  take  into  account  only  radiation  emission  (without  reabsorption)  this  temperature  will  be 
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approximately  600  K.  This  result  is  quite  natural,  because  losses  of  thermal  energy  by  heat  radiation  are 
rather  great.  It  is  easy  to  estimate  the  role  of  heat  radiation  reabsorption  comparing  temperature 
distributions  at  identical  time  moments  for  different  models  of  radiation  heat  transfer:  the  significant  part 
of  the  radiant  energy  does  not  leave  the  heated  volume,  reducing  total  losses  of  the  fireball  energy.  It  is 
obviously  that  the  total  energy  balance  makes  significant  effect  on  all  parameters  of  the  fireball. 

Fields  of  temperature  and  species  concentrations  presented  above  were  used  for  calculation  of  spectral 
signature  of  the  fireball  at  normal  direction  to  its  axis  of  symmetry.  Figures  6.60,  6.61  show  spectral 
signatures  in  W/(pm  sr)  at  consecutive  time  moments.  The  Monte-Carlo  simulation  technique  was  used  in 
180  spectral  bands  covered  spectral  range  AQ  =  1 000^-  9000  cm-1. 


Wavenumber,  1/cm 


(a) 


Wavenumber,  1/cm 


(b) 


Figure  6.60:  Spectral  signature  (W/(pmsr))  of  the  fireball  at  t  =  2  s  and  5  s 


Wavenumber,  1/cm 


(a) 


Figure  6.61:  Spectral  signature  (W/(pmsr))  of  the  fireball  at  t  =  11  s  and  15  s 


Homogeneous  spectral  grid  was  used  here.  One  million  photon  groups  were  modeled  in  the  each  spectral 
group.Analysis  of  the  data  presented  allows  conclude  that  time  dependence  of  spectral  signature  adequate 
reflects  facts  of  increasing  of  the  fireball  size  and  cooling  processes.  One  can  see  that  the  spectral 
signature  in  short-wave  part  of  investigated  spectral  region  (at  wave-number  larger  than  6000  cm1)  drops 
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more  quickly  than  in  long-wave  part  due  to  fast  cooling.  Spectral  emissivity  in  the  long-wave  part  (at 
wave-number  — >  1000  cm-1)  corresponds  to  relative  low-temperature  regions  of  the  fireball.  Therefore, 
practically  stable  radiation  in  this  spectral  region  is  clarified  by  the  following:  not  very  high  rate  of 
cooling  of  the  fireball  in  temperature  region  ~  400  4-  600  K  is  compensated  by  increasing  of  its  size. 

So,  presented  in  [51]  radiative  gas  dynamic  model  of  fireballs  generated  at  rocket  explosions  at  active  part 
of  trajectory  allows  predict  time-dependence  of  spectral  signature.  It  is  shown  that  to  predict  adequate 
spectral  signature  of  fireball  it  is  necessary  to  take  into  account  not  only  gas  dynamic  and  burning 
processes,  but  also  radiation  heat  transfer  processes.  Presented  numerical  simulation  results  show  that 
typical  time  for  registration  of  a  rocket  explosion  by  spectral  signature  is  more  than  ~30  s. 
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